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''        Major  Department:  Aerospace  Engineering,  Mechanics, 

and  Engineering  Science 

The  problem  of  time-dependent  enhanced  heat  transfer 

in  an  incompressible  viscous  laminar  fluid  subjected  to 

sinusoidal   oscillations   in   circular   pipes   which   are 

connected  to  a  hot  reservoir  at  one  end  and  a  cold  reservoir 

'        at  the  other  end  has  been  examined  numerically  in  detail. 

;.        Three  models  were  designed  for  the  investigation  of  such  an 

■^       enhanced  thermal  pumping  process  and  a  computer  code  (ETP) 

was  developed  to  implement  all  the  numerical  calculations. 

To  increase  the  understanding  of  the  mechanism  of 

;''       thermal   pumping,   the   periodic   velocity   profiles   and 

Lagrangian  displacements  as  well  as  tidal  displacements  at 

various  Wormersley  numbers  (from  a.     =  0.1  to  1000)  were 

studied.    Some  transient  problems  of  enhanced  axial  heat 

transfer  in  oscillating  pipe  flow  such  as  the  periodic  final 

temperature  build-up  process  in  oscillating  pipe  flow  were 

also  examined.   The  time-dependent  temperature  distribution 

xii 


in  the  different  models  was  numerically  studied  in  detail. 
The  enhanced  axial  heat  flux  magnitude  versus  different 
tidal  displacements  with  water  and  mercury  as  the  working 
fluids  bounded  by  pipe  walls  of  different  material  were 
observed   and   the   quadratic   coefficients   found.     The 
influence  of  the  variation  of  water  properties  on  the 
enhanced  axial  heat  flux  was  numerically  examined  and  the 
results  show  that  the  enhanced  axial  heat  flux  can  vary 
about  150  percent  even  within  the  temperature  range  from  CC 
to  100° C.   The  effects  of  wall  thickness  and  pipe  diameter 
in  enhanced  thermal  pumping  were  also  studied  and  the 
optimum  wall  thickness  was  found  to  be  about  20  percent  of 
the  pipe  radius  in  the  water-glass  combination.   The  tuning 
effect  in  the  water-steel  and  the  mercury-steel  cases  was 
examined  and  the  results  show  good  agreement  with  analytic 
predictions.   A  comparison  of  the  enhanced  axial  heat  flux 
with  the  axial  heat  flow  due  to  heat  conduction  at  various 
tidal  displacement  and  Wormersley  numbers  shows  that  the 
latter  is  quite  small  and  negligible  provided  the  tuning 
condition  is  satisfied. 

This  study  has  shown  that  the  enhanced  thermal  pump 
is  indeed  a  very  effective  tool  for  those  problems  where 
large  amounts  of  heat  must  be  transported  without  an 
accompanying  convective  mass  exchange.  The  investigation 
also  indicates  that  turbulent  flow  in  the  reservoirs  is 
preferable  to  laminar  conditions  and  should  receive  more 
attention  in  future  studies. 
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CHAPTER  I 
INTRODUCTION 


Enhanced  heat  transport  in  a  viscous  laminar  fluid 
subjected  to  sinusoidal  oscillations  in  a  very  long  pipe 
which  connects  a  hot  fluid  reservoir  at  one  end  and  a  cold 
fluid  reservoir  at  the  other  end  (Fig.  2-1)  has  been 
recognized  and  studied  recently  by  Kurzweg  [15,  16,  17,  22]. 
The  results  obtained  show  that  with  this  oscillatory  pipe 
flow  the  heat  transferred  axially  from  the  hot  end  to  the 
cold  end  can  be  orders  of  magnitude  larger  than  that 
obtained  by  pure  molecular  conduction  in  the  absence  of 
oscillations.  In  addition,  the  more  important  thing  of 
interest  is  that  this  heat  transfer  process  involves  no  net 
convective  mass  transport.  Major  assumptions  made  in  the 
above  cited  studies  on  enhanced  heat  diffusion  are  that  a 
constant  time-averaged  non-zero  axial  temperature  gradient 
is  always  present  in  the  oscillating  flow  and  that  the  axial 
molecular  conduction  along  the  wall  and  in  the  oscillating 
fluid  is  negligible. 

Discovery  of  this  enhanced  heat  transport  phenomenon 
was  made  possible  by  earlier  studies  on  axial  dispersion  of 
contaminants  within  steady  laminar  flows  through  capillary 
tubes  by  Taylor  in  1953  [32],  and  Aris  in  1956  [4].   These 
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earlier  studies  show  that  when  a  small  quantity  of  a 
contaminant  is  introduced  into  a  circular  pipe,  the 
dispersion  of  the  resultant  contaminant  cloud  is  greatly 
enhanced  by  the  flow  of  the  fluid.  Bowden's  1965  results 
^^y  show  that  similar  dispersion  effects  occur  in  oscillatory 

'    -V. . 

-     flow  [8].   This  enhanced  axial  dispersion  of  contaminants  in 
the  presence  of  laminar  oscillatory  flow  within  capillary 

p.|>  tubes  was  studied  in  1975  by  Chatwin  who  suggested  that  the 
assumption  of  constant  time-averaged  axial  contaminant 
gradient  can  be  made  [10].  Recently,  Bohn  et  al.  extended 
this  work  to  the  study  of  gas  component  transfer  in  binary 
gas  mixtures  when  these  are  confined  to  single  tubes  and  a 
sinusoidal  pressure  variation  is  applied  [7],  Further 
studies  in  1983  by  Watson  [38]  show  that  the  effective 
diffusion  of  contaminants  is  proportional  to  the  square  of 
the  tidal  displacement.  This  has  been  experimentally 
verified  by  Joshi  et  al.  [13]  and  by  Jaeger  [12],  both  in 
1983,  and  most  recently  by  Kurzweg  and  Jaeger  [19],  in  1987. 

*■  All  these  results  show  that  the  contaminant  would  spread 
axially  in  both  steady  and  oscillatory  laminar  pipe  flow  at 

»        rates  as  much  as  five  orders  of  magnitude  higher  than  in  the 

•'-  -4.      absence  of  fluid  motion. 

The  first  significant  research  work  extending  these 
enhanced  axially  diffusion  studies  to  the  heat  transfer 
problem  in  oscillatory  flow  within  very  slender  pipes  or 
flat  plate  channels  is  due  to  Kurzweg  [15,  16,  22].    In 
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early  1983  Kurzweg  suggested  that  a  similar  dispersion 
process  should  occur  in  the  heat  transfer  area  because  of 
the  similarity  in  both  the  governing  diffusion  and  heat 
transfer  equations  [16],  and  the  first  preliminary  theory 
was  formulated  in  1985  [15,  16],  in  which,  referring  to 
Chatwin's  idea,  a  time  averaged  constant  axial  temperature 
gradient  assumption  was  used.  The  instantaneous  temperature 
distribution  was  taken  to  be  of  the  form  [15,  22]. 

T  =  7[x  +  Ri  g(r)e^'"^]  (1_1) 

where  7  =  aT/ax  is  the  time-averaged  axial  temperature 
gradient,  R^  is  the  tube  radius,  x  is  the  axial  distance 
along  the  capillaries  (with  x  =  L/2  and  x  =  -L/2  denoting 
the  ends)  ,  L  is  the  pipe  length  under  consideration,  w  is 
the  oscillating  frequency  of  fluid,  g(j-)  is  a  radial 
temperature  distribution,  and  f  =  r/Ri  is  the  dimensionless 
radial  distance.  The  theoretical  analysis  shows  that  under 
certain  conditions  such  enhanced  axial  heat  diffusivity  can 
indeed  be  significantly  larger  than  the  axial  molecular 
conduction  [15,  16]  and  this  has  been  verified  by  some 
experimental  measurements  by  Kurzweg  and  Zhao  [22]. 

In  order  to  better  understand  the  physical  mechanisms 
of  this  interesting  and  potentially  useful  heat  transfer 
technique,  we  shall  examine  in  greater  detail  the  thermal 
pumping  model  shown  in  Fig.  2-I.  it  is  assumed  that  a 
bundle  of  very  thin  and  long  tubes  is  connected  to  a 


4 
reservoir  which  supplies  unlimited  hot  liquid  at  one  end  to 
a  second  reservoir  which  supplies  unlimited  cold  liquid  at 
the  other  end.    The  liquid  in  the  pipes  is  oscillating 
axially  with  an  amplitude  such  that  none  of  the  liquid  which 
is  originally  in  the  middle  portion  of  pipes  ever  runs  into 
either  reservoir.   That  is,  there  is  no  net  convective  mass 
exchange  between  two  reservoirs.   The  largest  axial  fluid 
dimensionless  displacement  (when  cross-section  averaged)  is 
referred  to  as  the  nondimensional  tidal  displacement  and  is 
denoted  by   "AX"   (it   should  not  be  confused  with  the 
dimensional  tidal  displacement  "ax"  frequently  used  in  the 
present  study).   At  time  t  =  0,  the  fluid  within  the  pipes 
is  set  into  axial  oscillations  at  angular  frequency  u    and 
tidal  displacement  Ax.    After  a  short  transient,   this 
oscillatory  motion  will  lead  to  very  large  axial  heat  flows 
which  can  be  readily  made  to  exceed  those  possible  with  heat 
pipes. 

Before  exploring  the  mechanisms  of  this  enhanced  heat 
transport,  it  is  necessary  to  introduce  some  new  concepts 
which  are  commonly  used  in  the  study  of  this  type  of 
oscillatory  motion. 

As  is  well  known  [36],  for  high  frequency  viscous 
laminar  axial  oscillations  within  fluid  flow  along  rigid 
pipes,  the  non-slip  boundary  condition  creates  a  very  thin 
Stokes'  viscous  boundary  layer  of  thickness 

S   =  f2U7 


^t-  °  -  V^^/<^  (1-2) 
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where  u  is  the  fluid  kinematic  viscosity.  For  room 
temperature  water  at  a  frequency  of  10  Hz,  this  viscous 
boundary  layer  is  approximately  1.7.10"2  cm.  The 
corresponding  thermal  boundary  layer  thickness  is  about 

5th  =  5//Pr  (1_3) 

where  Pr  is  the  Prandtl  number.   Note  that  both  8    and  5th 
decrease   in   thickness   as   the   oscillating   frequency  w 
increases. 

In  the  theoretical  analysis,  it  is  always  assumed 
that  a  fully  developed  velocity  profile  of  the  oscillating 
flow  exists  within  the  pipes.  At  high  frequency  w,  this 
flow  consists  essentially  of  a  slug  flow  over  most  of  the 
fluid  core  bounded  by  a  thin  boundary  layer  of  width  8  as 
discussed  by  Uchida  [36].  Neglecting  end  effects,  the  fully 
developed  oscillating  laminar  velocity  profile  in  pipes,  due 
to  a  periodic  axial  pressure  gradient,  is  found  to  be  [19] 


U(r,t)  =  Uof(r)e^'"'^  (1-4) 


where  Uq  is  a  representative  velocity,  $•  =  r/R^  is  again  the 
dimensionless  radial  distance,  f($-)  the  velocity  shape 
function,  and  w  the  angular  velocity  of  the  oscillatory 
flow.   The  explicit  form  of  f($-)  is 


:^:.': 


'  ?-<> 


f(0  =  -~2r-   [   1  -  ^^H"-"P   1  (1-5) 
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where  A  =  R^^  |ap/ax|/pi/Uo  is  the  nondimensional  pressure 
gradient  maximum  acting  along  the  capillaries,  a  =  R^^/u/u  is 
the  Wormersley  number  measuring  the  ratio  of  inertia  to 
viscous  forces,  i/  is  the  fluid  kinematic  viscosity  of  the 
carrier  liquid,  and  p  is  the  fluid  density.  This  velocity 
profile  will  reduce  to  the  familiar  Poiseuille  parabolic 
shape  as  the  angular  frequency  w  becomes  small,  while  at 
moderate  frequency,  f($-)  has  the  shape  demonstrated  by 
Uchida  [36]. 

Another  new  term  commonly  used  when  dealing  with 
oscillating  flow  is  the  cross-stream  averaged  dimensionless 
tidal  displacement  AX  which  can  be  mathematically  defined  as 


4Uo   I  rl 


^x  =  — z—      Lrf(r)dr 


and  on  integration,  yields  [18] 


(1-6) 


|ap/ax| 

AX  = 


1/2  pw^ 


1  +  -^-   F(a) 


(1-7) 


Where  the  complex  function  F(a)  has  the  form 


F(a)  =  i   ^n'(/^) 


with   the   prime   denoting   differentiation.     Using   the 
definition  of  Kelvin  function:   Jo(y^i^)  =  ber  a  +  i  bei  a 
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[1],  the  complex  function  F(q)   (Fig.  l-i)  can  be  further 
written  as: 

F(«)  =  F,(«, .  i  Fi«.,  =  i  [  III': :  i  a':  ]      d-s) 

This  dimensionless  tidal  displacement  is  related  to  the 
maximum  of  the  periodic  pressure  gradient  via  [18] 


f 

|ap/ax| 

AX  =  T-7T 2- 

1/2  pu>^ 


a'^  a  «  1 
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/2  (1-9) 

1 a  »  1 


Apparently,   for   any   fixed   tidal   displacement   ax   and 
oscillating   frequency   w,   the   axial   pressure   gradient 
required  is  proportional  to  the  inverse  square  of  the 
Wormersley   number   when  a       is   small;   however,   it   is 
independent  of  a    when  a    is  very  large.   This  also  implies 
that   the   exciting   axial   pressure   gradient   |ap/ax|   is 
approximately  proportional  to  the  fluid  kinematic  viscosity 
1/  and  inversely  proportional  to  the  square  of  the  pipe 
radius  R^  when  Wormersley  number  is  very  small  (it  happens 
only  at  low  w,  small  R^,  and  large  u,     for  example,  oil), 
while  it  is  almost  independent  of  the  fluid  kinematic 
:       viscosity  and  the  pipe  radius  when  the  Wormersley  is  very 

V'-'tr 

||.;:     large  (it  happens  only  at  very  large  R^,  high  w,    small  u, 
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for  instance,  a  liquid  metal) .  Note,  if  the  tidal 
displacement  AX  is  fixed,  while  allowing  the  oscillating 
frequency  to  change,  the  axial  pressure  gradient  |ap/3x|  can 
become  very  large  when  the  oscillating  frequency  becomes 
large.  This  is  mainly  due  to  inertial  effects  and  not  so 
much  due  to  viscous  drag  forces  which  dominate  the 
oscillatory  flow  at  small  Wormersley  number. 

With  the  definition  of  the  above  quantities,  we  are 
now  in  a  position  to  explore  the  details  of  the  enhanced 
axial  heat  transfer  in  oscillatory  flow  within  pipes.   It  is 
assumed  that  a  constant  temperature  gradient  exists  along 
the  pipe  in  the  axial  direction  and  that  a  very  large  time- 
dependent   radial   temperature   gradient   variation   is 
superimposed.   When  the  fluid  moves  towards  the  cooler  side 
(we  term  this  the  positive  stroke) ,  the  hotter  fluid  within 
the  pipe  core  which  is  initially  brought  into  the  pipe  from 
the  hot  reservoir  produces  a  large  radial  heat  flow  via 
conduction  to  the  cooler  portions  of  the  fluid  within  the 
Stokes-  boundary  layer  and  to  the  cooler  pipe  wall;  while 
during  the  negative  (or  reverse)  stroke,  i.e.,   when  the 
fluid  moves  towards  the  hotter  side,  the  higher  temperature 
in  the  boundary  layer  and  the  pipe  wall  conducts  the  heat 
back  into  the  cooler  fluid  core.  This  coupled  radial  heat 
conduction  with  an  axial  convective  transport  leads  to  an 
enhanced  axial  heat  flux  along  the  entire  length  of  the 
pipes. 
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Further,  from  the  system  point  of  view,  the  heated 
fluid  near  the  cold  reservoir  will  eventually  be  ejected 
into  the  cold  reservoir  and  mixed  there  with  the  lower 
temperature  liquid.   Contrarily,  near  the  hot  reservoir  side 
the  fluid  within  the  pipes  which  has  been  cooled  during  each 
positive  stroke  is  pushed  out  into  the  hot  reservoir  and 
mixes  with  the  higher  temperature  liquid.   This  process  of 
thermal  pumping  is  what  leads  to  a  time-averaged  heat  flow 
from  the  hot  reservoir  into  the  cold  reservoir.   It  differs 
essentially  from  the  working  principle  of  a  normal  pump. 
For  a  normal  pump  one  can  draw  the  analogy  with  transport  of 
a  one-way  vehicle  which  transports  passengers  as  well  as  the 
carrier  from  one  point  to  another.   For  the  thermal  pump  one 
can   draw   the   analogy   with   a   two-way   busline   which 
periodically  loads  and  unloads  the  passengers  (heat)  from 
the  hot  reservoir  to  the  cold  reservoir,   heat  can  be 
continually  transferred,   and  the  carrier,   in  the  time- 
averaged   sense,   does   not   move.     This   property   is 
particularly  important  for  those  systems  where  a  large 
amount  of  heat  transfer  is  needed  while  the  working  fluid  is 
required  to  remain  in  the  system  (as  in  nuclear  reactors) . 
Note  that  the  axial  heat  conduction,  in  general,  is  assumed 
to  be  very  small  in  this  thermal  pumping  process  compared 
with  the  enhanced  axial  heat  flow  [16]. 

Apparently,   the  heat   transport   rate   in   thermal 
pumping  is  governed  by  both  the  thermal  properties  of  the 
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working  medium  and  pipe  wall  and  the  characteristics  of  the 
oscillatory  motion  of  the  fluid.   The  enhanced  axial  heat 
flow  does  increase  with  increasing  oscillating  frequency  as 
this   thins   out   the   boundary   layer   and   leads   to   an 
accompanying  increase  in  the  radial  temperature  gradient. 
This   observation   holds   only   as   long   as   the   thermal 
properties  of  the  liquid  and  the  wall  are  compatible.   If 
the  molecular   conduction   of   the   fluid   in  the  radial 
direction   is   very   small,   then   even   high   frequency 
oscillatory  motion  will  not  produce  a  large  increase  in  the 
rate  of  enhanced  axial  heat  flow.   This  is  because  such  a 
system  fails  to  supply  the  "passengers"  enough  time  to  be 
loaded  onto  the  "bus"  and  to  be  unloaded  from  the  "bus". 
Obviously,  the  system  just  wastes  energy.    On  the  other 
hand,  if  the  molecular  conduction  of  the  fluid  in  the  radial 
direction   is   very   large,   but   the   frequency   of   the 
oscillatory  motion  is  low,  once  again  one  can  not  expect 
that  there  will  be  an  efficient  enhanced  axial  heat  transfer 
between  two  reservoirs  because  the  "bus"  is  now  moving  too 
slowly. 

The  above  observation  can  be  confirmed  by  an  analysis 
of  the  performance  of  a  water-glass  combination  (i.e.,  water 
is  the  working  fluid  medium  within  a  glass  tube)  and  of  a 
mercury-steel  combination.  For  the  former,  it  is  necessary 
to  employ  rather  small  diameter  tubes  and  low  frequency 
oscillations  with  large  tidal  displacement,   for  it  has 
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relatively  poor  heat  diffusive  properties  as  compared  with 
the  mercury-steel  case.    This  small  tube  diameter-lower 
oscillation  frequency  set-up  is  necessary  in  order  to  ensure 
that  there  is  sufficient  time  to  transfer  the  excessive  heat 
content  of  the  bulk  water  core  to  the  tube  wall  during  the 
positive  stroke  of  each  period  and  to  permit  the  transfer  of 
excess  heat  from  wall  to  the  cooler  core  fluid  during  the 
negative  stroke.   Otherwise,  the  water  would  just  carry  a 
portion  of  its  heat  content  back  and  forth  in  the  pipe  and 
the  condition  for  achieving  optimum  enhanced  axial  heat 
transfer  could  not  be  met.   For  the  mercury-steel  case,  one 
can   chose   relatively   large   pipe   diameters   and  higher 
oscillation   frequencies  with  smaller  tidal   displacement 
because  of  the  higher  thermal  conductivity.   This  assures 
that  only  very  short  times  are  needed  to  exchange  the  heat 
between  the  core  of  the  fluid  and  the  wall.   it  should  be 
pointed  out  that  the  suggestion  of  using  smaller  tidal 
displacement  is  purely  due  to  the  mechanical  considerations 
and  that  one  always  tries  to  keep  ax  as  large  as  possible  in 
order  to  produce  large  axial  heat  flows. 

From  the  above  discussion,  it  can  be  concluded  that 
the  process  of  enhanced  heat  transfer  via  oscillatory 
pumping  requires  a  precise  tuning  of  parameters  governing 
the  enhanced  heat  transport.  Indeed  there  is  an  expected 
"tuning  effect"  as  discussed  in  references  [16,  20,  21]. 
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The  tuning  effect  is  a  very  important  concept  in  the 

study  of  the  presently  considered  heat  transfer  process.   It 

shows  that  there  will  be  an  optimum  combination  between 

thermal  properties  of  the  working  medium  and  wall  and  the 

characteristics  of  the  oscillatory  motion.   The  qualitative 

aspects  of  the  tuning  effect  have  been  observed  earlier  for 

both  the  case  of  a  flat  channel  and  that  of  the  cylindrical 

pipe  [15,  16].   From  Fig.  4-38,  one  can  see  that  an  optimum 

for  axial  heat  transfer  occurs  only  at  or  near  the  tuning 

point  which  depends  on  the  oscillating  frequency  and  the 

thermodynamic  properties  of  the  fluid  and  wall.   As  has  been 

pointed  out  by  Kurzweg  [20],  in  order  to  obtain  the  optimum 

enhanced  heat  transfer  one  has  to  carefully  select  suitable 

values  for  the  pipe  size,  the  material  of  pipe  and  the 

working  medium  as  well  as  the  manner  of  oscillatory  motion. 

The  nondimensionalized  enhanced  heat  diffusivity  is 
defined  as 


A  = 


Ke 


wAX^  (1-10) 


Where  Kg  =  4>/ipc  is  the  coefficient  of  enhanced  heat 
diffusivity,  <!>  is  the  axial  heat  flux,  7  is  the  time- 
averaged  axial  temperature  gradient,  p  is  the  density  of  the 
fluid,  and  c  is  the  specific  heat.  One  can  show  that  this 
nondimensional  enhanced  heat  diffusivity  is  a  function  of 
both  the  Wormersley  number  and  the  Prandtl  number  [21]  and 
hence  that  the  dimensional  axial  heat  diffusivity  Kg  is  a 


-J 


14 

function  of  the  tube  radius  Ri,  the  oscillating  frequency  w, 
the  kinematic  viscosity  u,  and  the  square  of  the 
dimensionless  tidal  displacement  AX.  This  can  be  explained 
from  the  fact  that  the  radial  heat  flow  is  proportional  to 
the  product  of  the  representative  radial  temperature 
gradient  7AX/5th  and  the  surface  area  per  unit  depth  of  ttAX 
available  for  cross-stream  heat  transport. 

The  use  of  large  tidal  displacement  is  always 
beneficial  in  the  enhanced  axial  heat  transfer  within 
oscillating  pipe  flow.  However,  in  order  to  avoid  the 
direct  convective  net  mass  exchange  between  the  two 
reservoirs,  the  tidal  displacement  must  be  limited  to  less 
than  about  one  half  of  the  pipe  length. 

As  has  already  been  predicted  by  theoretical  studies 
and  will  be  confirmed  by  the  present  numerical  simulations, 
the  axial  heat  transfer  will  be  further  enhanced  if  the 
rigid  surface  (part  of  the  rigid  wall  with  finite  thickness) 
has  a  non-zero  thermal  diffusivity  and  hence  heat  storage 
capability. 

Note  that  the  existing  considerations  are  restricted 
to  laminar  flow.  Turbulent  flow  conditions  can  occur  in 
oscillating  pipe  flow  at  higher  values  of  o.AxVi'  [26,  27] 
and  apparently  would  destroy  the  assumptions  of  the  current 
analytic  model  of  the  thermal  pumping  process.  Fortunately, 
the  condition  for  optimum  enhanced  heat  transfer  in  such 
oscillating  pipe  flow  obtained  at  the  tuning  point  requires 
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very  slender  pipes,  such  that  the  Reynolds  number  is  usually 
small  enough  so  that  the  oscillating  motion  falls  within  the 
laminar  range  [22]. 

The  theoretical  aspects  of  the  oscillatory  enhanced 
axial  heat  transfer  process  have  been  developed  much  further 
than  its  experimental  and  numerical  counterparts.  The 
theoretical  predictions  are  quite  limited  and  consider  only 
cases  under  certain  simplifying  assumptions  [8].  Numerical 
work  is  necessary  in  order  to  not  only  to  examine  the 
correctness  of  the  theoretical  analysis  but  also  to  further 
the  development  of  advanced  enhanced  thermal  pumping 
devices.  Numerical  studies  are  not  only  fast,  economical 
and  accurate,  but  they  also  offer  a  handy  way  to  access 
complex  geometries  which  can  not  be  handled  analytically. 

It  is  the  purpose  of  this  study  to  extend  the 
analytic  work  on  thermal  pumping  by  a  detailed  numerical 
study.  We  intend  first  to  examine  some  transient  problems 
of  axial  heat  transfer  in  oscillatory  pipe  flow,  such  as  the 
development  of  the  velocity  profile  at  various  Wormersley 
numbers  in  contrast  with  those  of  reference  [36],  where  only 
several  special  cases  with  intermediate  Wormersley  number  a 
were  discussed,  and  to  examine  numerically  the  relationship 
between  the  tidal  displacement  and  the  required 
corresponding  strength  of  the  periodic  pressure  gradient  as 
a  function  of  Wormersley  number  a    and  of  tidal  displacement 

AX. 
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Next,  we  examine  the  build-up  process  of  the 
temperature  distribution  in  a  pipe  which  connects  a  hot 
reservoir  at  one  end  to  a  cold  reservoir  at  the  other  end 
and  see  whether  there  actually  exists  a  constant  time- 
averaged  temperature  gradient  along  the  pipe  axis  when  the 
final  periodic  state  is  eventually  reached.  Note  that  a 
time-averaged  linear  temperature  distribution  along  the 
axial  direction  is  an  essential  assumption  in  the  existing 
theoretical  studies. 

The  third  part  of  this  investigation  which  forms  the 
main  effort,  is  a  computer-aided  numerical  simulation  of  the 
thermal  pumping  technigue,  including  an  investigation  of  the 
variation  of  the  enhanced  axial  heat  flux  versus  the  tidal 
displacement,  the  variation  of  enhanced  axial  heat  flux 
versus  different  Wormersley  numbers,  and  a  study  of  the 
variation  of  heat  flux  versus  different  Prandtl  numbers.  It 
also  includes  a  study  of  the  influence  of  wall  thickness  and 
pipe  diameter  as  well  as  the  change  of  the  fluid  properties 
on  such  an  enhanced  axial  heat  flux  and  an  examination  of 
the  tuning  effect  in  the  conducting  wall  case.   Further,  we 
compare  axial  heat  transfer  in  oscillating  flow  with  that  in 
the  steady  flow,  and  compare  quantitatively  this  enhanced 
heat  flux  with  that  induced  by  the  pure  axial  molecular 
conduction  under  various  Wormersley  numbers  and  different 
fluid-solid  wall  combinations.  Finally,  the  author  examines 
the  phase  lag  phenomenon  in  oscillatory  flow  in  order  to 
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further  reveal  the  mechanism  whereby  this  new  heat  transfer 
technique  functions. 

The  computational  work  presented  here  was  completed 
on  the  Vax-11/750  electronic  computer  in  the  Department  of 
Aerospace  Engineering,  Mechanics,  and  Engineering  Science, 
and  on  Vax-11/780   in  the  Center  for  Instructional  and 
Research   Computing   Activities,   University   of   Florida, 
Gainesville,   Florida.    The  numerical  approach  used,   for 
solving  the  presently  considered  heat  transfer  problem, 
employed   a   second-order   Crank-Nicolson   scheme   and   an 
Alternating  Direction  Implicit  method   (ADI)   with  Thomas 
algorithm.    A  Fortran  computer  code  named  ETP  (Enhanced 
Thermal   Pumping)   was   developed   to   implement   all   the 
calculations. 


CHAPTER  II 
FORMULATION  OF  THE  PROBLEM 

A  single  tube  with  inner  radius  R^,  outer  radius  R2 
and  length  L  (such  that  L  »  R^  and  R2)  connects  a  large 
reservoir  of  hot  fluid  at  temperature  T  =  Th  at  one  end,  to 
another   large   reservoir   of   the   same   fluid   at   cold 
temperature  T  =  Tc  at  the  other  end  (Figs.  2-1  and  2-2). 
The  tidal  displacement  ax  and  the  oscillating  frequency  w 
are  adjusted  to  assure  that  the  turbulent  flow  does  not 
occur.  The  tube  is  oriented  in  such  a  manner  so  that  the 
effect  of  gravity  on  the  oscillatory  motion  of  the  fluid  in 
the  pipe  is  negligible.   Water  and  mercury  are  employed  as 
the  working  mediums.   They  are  understood  to  be  Newtonian 
and   incompressible   fluids.     The  variations  of  thermal 
properties  of  the  fluid  with  temperature  during  the  heat 
transport  process  are  assumed  to  be  negligible.   With  these 
assumptions  the  problem  of  enhanced  time-dependent  heat 
transfer  induced  by  a  simple  harmonic  oscillatory  laminar 
fluid  motion  in  a  very  long  circular  pipe  can  then  readily 
be  formulated. 

Governing  Equations 
The  use  of  a  very  slender  capillary  circular  tube 
with  constant  cross-sectional   area   and  neglecting  end 
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effects  insures  a  laminar  axisymmetrical  one-dimensional 
time-dependent  motion.  It  is  convenient  to  employ 
cylindrical  coordinates  for  this  problem  and  we  denote  the 
coordinate  in  the  axial  direction  of  the  tube  by  x,  and  the 
radial  direction  by  r.  The  axial  velocity  can  be  taken  to 
be  independent  of  x  and,  in  order  to  satisfy  the 
requirements  of  continuity  the  other  velocity  components 
must  vanish.  We  shall  further  assume  that  the  pressure 
gradient  induced  by  moving  the  piston  (Fig.  2-1)  is  harmonic 
and  has  the  form  [29] 

1    8p 

"  ~r  "air  =  ^  """^  ""^  (2-1) 

where  A  =  l/p\dp/dx\  is  a  constant  which  measures  the 
maximum  pressure  gradient  existing  along  the  x-axis. 
Clearly  this  equation  implies  that  we  are  now  dealing  with  a 
time-dependent  sinusoidal  pressure  gradient  which  is 
constant  over  the  pipe  cross-section  at  any  instant  and  the 
pressure  varies  linearly  along  the  x-axis.  The  simplified 
Navier-Stokes  equation  for  this  problem  is  [36] 


au  u     a 


d 


du 


^  ^  u  r  (J  \J  -I 

rr~  =  Acos(wt)  +  - —  r(z — ) 

t        ^   '    r  ar   L  ^dr  ' j 

0  <  r  <  Ri 


(2-2) 


where  U(r,t)   is  the  time  and  radially  dependent  axial 
velocity  component. 
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The  corresponding  temperature  T(x,r,t)  of  the  fluid 
within  the  pipe  is  governed  by  the  heat  conduction  equation 
[15] 


3T      aT     p  1   a     aT     s^t 


r   ar  ^^  ar  ^  "*"  dx'^     J 

0  <  r  <  Ri 

(2-3) 

where  R^  is  the  inner  radius  of  the  tube  and  ^f  the  fluid 
thermal  diffusivity  which  is  related  the  thermal 
conductivity  kf  by 

kf 


Kf 


pC 


Here  p  is  the  fluid  density,  and  c  is  the  specific  heat. 
Note  that  the  viscous  heating  term  has  been  neglected  in 
equation  (2-3)  since  it  is  very  small  for  most  experimental 
conditions;  this  is  justified  provided  one  does  not  deal 
with  very  high  Prandtl  number  fluids  such  as  oils.  The 
temperature  in  the  wall  can  be  determined  from  the  solution 
of 

ST  r  "'"    ^        ^"^        ^^'^  1 

Tt'  =  ""^[V  a^('^~a^)  +  T^ 

Rl  <  r  <  R2 

(2-4) 
where  /c^  is  the  thermal  diffusivity  of  the  conducting  wall 
in  Rl  <  r  <  R2  and  is  defined  by 
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K-^     - 


where  k^^,  p^,  and  c^,  are  the  wall  conductivity,  density, 
and  the  specific  heat,  respectively.  By  introducing  the 
following  non-dimensional  terms,  equation  (2-2) ,  (2-3)  and 
(2-4)  can  be  treated  more  easily. 


t*  = 


r*  = 


1/w 

r 


X   = 


u*  = 


s 

u 


A/w 


m*   ^ 


Th  -  Tc 


Where  6  =  JTvJZ)  is  again  the  fluid  viscous  boundary  layer 
thickness.  The  dimensionless  governing  system  of  equations 
can  then  be  written  as 


au' 


at' 


aT' 


at 


and 


aT' 


at' 


=  cos(t*)  + 


IT-   =    CU' 


2S 


a^u* 


au' 


ar 


w^~  + 


ar' 


aT' 


ax' 


1  r  a2' 
!Pr  [~dr 


a^T* 
W2~  + 


2ip* 


aT' 


0  <  r  <  Ri 


(2-5) 


aT' 


a^T* 


Ic^     dr 


w  + 


0  <  r  <  R^ 


(2-6) 


7^ 


ar' 


a  2rn* 

Rl  <  r  <  R2 


(2-7) 


23 


where 


Pr  =  

K, 


s  = 


AC 


w 


and 


-A 


This  nondimensionalization  has  some  advantages  in  the 
computing  process  to  be  carried  out  below.  The 
dimensionless  velocity  and  its  distribution  over  the  cross- 
section  found  from  the  momentum  equation  (2-5)  are  expected 
to  be  universal  for  any  Wormersley  number  a  =  R^Jup/n  and 
any  associated  quantities,  such  as  the  tidal  displacement 
and  Lagrangian  displacement.  Its  final  periodic  form  is  of 
the  form  given  by  (1-4).  The  dimensionless  temperature  in 
the  pipe  is  only  related  to  the  Prantdl  number  and  the 
dimensionless  velocity,  while  that  in  the  wall  is  related  to 
the  ratio  of  wall  heat  diffusivity  to  the  kinematic 
viscosity,  as  seen  from  equation  (2-6)  and  (2-7) . 

The  governing  dimensionless  equations  (2-5) ,  (2-6) , 
and  (2-7)  are  a  set  of  second-order  parabolic  type  of 
partial  differential  equations  expressed  in  cylindrical 
coordinates.  To  solve  this  set  of  simultaneous  equations,  a 
corresponding  set  of  boundary  conditions  and  initial 
conditions  are  required.  Since  the  velocity  is  assumed  to 
be  a  function  of  r  and  t  only,  just  two  boundary  conditions 
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are  needed  for  solving  the  momentum  equation,  while  for  the 
temperature  T(x,r,t)  the  heat  conduction  equations  require 
two  boundary  conditions  in  both  the  r  and  x  directions  as 
well  as  compatibility  conditions  along  the  interface  between 
the  fluid  and  the  solid  wall.  It  should  be  pointed  out  that 
the  initial  conditions  are  less  important  than  the  boundary 
conditions  if  one  seeks  only  a  final  periodic  state. 

Boundary  Conditions 
The  boundary  conditions  for  the  velocity  are  the 
usual  ones  for  viscous  flow,  namely,  that  the  velocity  is 
zero  at  the  inner  surface  of  the  wall  (r*  =  R^*)  .  Also  by 
symmetry,  the  normal  derivative  of  velocity  at  the  axis  is 
zero.   That  is, 

U*(Ri*,  t*)  =  0  (2-8) 


and 


au*(0,  t*) 

JF"  =  °  (2-9) 

The  boundary  conditions  for  the  temperature  depend  on  the 
particular  model  investigated. 
Model  1 

In  this  model,  it  is  assumed  that  the  temperature  at 
each  end  of  the  pipe  is  equal  to  that  of  the  connecting 
fluid  filling  the  reservoirs  when  the  fluid  moves  into  the 
capillary  tube  (i.e.,  during  one  half  of  each  cycle).   The 
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end  boundary  temperature  at  x  =  0  and  x  =  L  will  take  the 
values  of  temperature  of  the  adjacent  pipe  fluid  elements. 
With  this  model,  we  tried  to  simulate  the  real  experimental 
observations  where  it  can  be  clearly  seen  that  a  fluid  jet 
exits  the  pipe  during  the  outstroke,  while  the  elements  of 
the  fluid  enter  the  tube  in  a  funnel  pattern  during  the 
instroke.  This  observation  shows  that  there  is  enough  time 
during  each  oscillating  cycle  to  allow  the  fluid  elements 
within  the  exiting  jet  to  fully  mix  with  the  fluid  particles 
which  originally  are  in  the  reservoir  before  the  next 
instroke,  provided  the  oscillating  frequency  is  not  too 
high.  The  end  temperature  boundary  conditions  are  here 
taken  as 


1* 
f 

T''(0,  r*,  t*)  = 


'"^  h  (when  fluid  enters  pipe) 


T*(L*,  r*,  t*)  = 


'^   adj         (when  fluid  leaves  pipe) 


(2-10) 

'"^^c  (when  fluid  enters  pipe) 

■^  adj        (when  fluid  leaves  pipe) 


■,* 


(2-11) 


where   T  c/  T  ^  are  the  nondimensional  temperatures  of  cold 
and  hot  reservoirs,   respectively.    While  T*adj   is  the 
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temperature  of  the  fluid  elements  which  are  adjacent  to  the 
corresponding  ends  at  a  particular  instant. 

The  thermal  boundary  conditions  along  the  outer 
surface  r*  =  R2  of  the  pipe  wall  are  taken  to  satisfy  the 
insulating  wall  condition,  and  along  the  axis  of  the  tube 
the  temperature  is  assumed  to  meet  the  symmetric  boundary 
requirement,  i.e.,  the  radial  derivatives  of  temperature 
along  axis  are  equal  to  zero.  We  thus  have 
^* 


dT" 


dr" 


r*  =  R2* 


=  0 


(2-12) 


and 


aT' 


ar' 


r*  =  0 


=  0 


(2-13) 


The  compatibility  conditions  along  the  interface 
between  the  fluid  and  the  solid  wall  are  that  the  radial 
heat  flux  and  temperature  are  constant  across  the  interface. 
That  is 


^f 


and 


aT' 


ar' 


aT' 


fluid 


-  kw  -^pr 


wall 


at  r  =  Ri 


(2-14) 


T   I        =  T 

'  fluid       wall 


at  r  =  Ri 


(2-15) 


Since  numerically  the  same  nodes  are  chosen  along  the 
interface,  condition  (2-15)  will  be  automatically  satisfied. 
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Model  2 

Here  it  is  assumed  that  a  heat  source  rim  of  width  2b 
is  directly  mounted  on  the  interface  between  the  solid  wall 
and  the  fluid  at  x  =  L/2,  while  two  small  cold  rim  sources 
of  width  b  each  are  mounted  at  x  =  0  and  x  =  L  (see  Fig.  2- 
3)  .  The  wall  thickness  is  assumed  to  be  zero.  This  model 
is  intended  to  simulate  the  enhanced  heat  transfer  process 
of  oscillatory  fluid  in  an  infinitely  long  pipe  which  is 
heated  and  cooled  by  the  alternative  evenly  distributed  heat 
and  cold  sources  along  the  outer  surface  of  a  very  thin  wall 
which  possesses  infinite  heat  conductivity.  Apparently, 
this  geometry  can  be  well  simulated  by  the  present  model 
with  periodic  boundary  conditions.  Mathematically  the 
boundary  conditions  are 

T*(x*,  Rl*,  t*)  =  T*h  X2*  <  X*  >  X3*      (2-14) 


0  <  X*  <  XI* 

T''(x'^,  Rl",  t")  =  T*c  or  (2-15) 

X4*  <  X*  <  L* 


3T*  XI*  <  X*  <  X2* 


dr 


^-  -  0  o^    ^  (2-16) 


X3*  <  X*  <  X4* 


where  XI  ,  X2*,  X3*,  X4*,  X5*  and  X6*  are  nondimensional 
coordinates  of  points  which  correspond  to  X  =  0,  b,  L/2-b, 
L/2+b,  L-b,  and  L,  respectively  in  Fig.  2-3.  The  periodic 
end  boundary  conditions  are  given  as 
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T*(0,  r*,  t*)  =  T*(L*,  r*,  t*)  (2-17) 


and 


dT*     I  dT* 


"^^  I  X*  =  0      "^ 


*     *  (2-18) 

X*  =  L* 


Model  3 

In  model  1,  the  boundary  conditions  for  the  heat 
equation  at  X  =  0  and  L  are  somewhat  artificial.    The 
temperature  distribution  in  the  vicinity  of  pipe  ends  may 
not  precisely  match  the  constant  temperature  end  conditions. 
This  is  particularly  true  in  the  large  tidal  displacement 
and/or   very   high   oscillating   frequency   cases   since   a 
temperature  drop  occurs  at  the  hot  end  and  a  rise  at  the 
cold  end  as  fluid  particles  exit  the  pipe.   This  leads  to  a 
discontinuity  of  temperature  which  will  lead  to  unexpected 
dispersion  errors  if  an  even  order  numerical  method  is  used 
[3,   24,   30].     To   improve  the  end  boundary  conditions 
presented  in  model  1,  the  following  extended  model  has  been 
considered.   It  is  assumed  that  an  extension  pipe  which  is 
some  5  times  the  length  of  the  central  pipe  is  attached  to 
the  original  pipe  used  in  model  1.   The  heat  and  cold  source 
are  assumed  to  be  located  on  the  outer  surface  of  the  wall 
in  the  extended  sections  of  the  pipe  as  well  as  at  both  ends 
of  the  tube  as  shown  in  Fig.  2-4.   This  model  is  used  to 
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investigate  the  situation  where  the  heat  contained  in  the 
jet  is  exchanged  by  pure  heat  conduction  with  the 
surrounding  fluid  elements  in  each  reservoir  without  any 
convective  mixing.  In  this  case,  the  boundary  conditions 
along  the  outer  surface  of  the  wall  are 


"^     -  '^  h    >  0<x*<  XI*  (2-19) 


T*  =  T*c  X2-  <  x"  <  X3«  (2-20) 

and 


dT* 


7* 


=    0  XI*  <  X*  <  X2 


ar^         "  A±   <  X   <  X2"  (2-21) 

where  XI  ,  X2*,  and  X3*  are  nondimensional  coordinates  of 
points  which  correspond  to  X  =  5L,  6L,  and  llL, 
respectively,  in  Fig.  2-4.   At  both  ends,  we  have 


T*(0,  r*,  t*)  =  T\  (2.22) 

and 


T*(L*,  r*,  t*)  =  T*c  (2-23) 


As  in  the  other  models,  the  symmetry  condition  along  the 
axis  of  pipe  requires  that 


dT* 


ar^   |r*=0   °  (2-24) 
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and  the  radial  flux  and  temperature  continuity  condition 
along  the  interface  are 


dT 


kf 


ar' 


and 


aT' 


fluid 


-  kw  -jp- 


wall 


(2-25) 


=  T' 


fluid 


wall 


(2-26) 


Here  again  kf  and  k^  are  the  same  as  described  in  the 
previous  models. 

Initial  Conditions 
The  initial  condition  chosen  depends  on  the  problem 
under  consideration.  If  one's  goal  is  to  investigate  the 
periodic  quasi-steady  solution  only,  the  initial  conditions 
chosen  should  be  as  close  to  the  quasi-steady  state  as 
possible  so  that  a  rapid  convergence  rate  at  low  CPU  time 
cost  is  achievable.  If  one  intends  to  study  the  transient 
process,  then  various  initial  conditions  should  be  supplied 
according  to  purpose  of  the  investigating  cases  selected. 
For  the  velocity  initial  condition  we  choose  for  all  our 
studies 


U*(x*,  r*,  0)  =  0 


(2-27) 


For  the  temperature,  if  the  purpose  of  investigation  is  to 
examine  the  build-up  process  of  the  periodic  quasi-steady 
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state  in  the  thermal  field,  the  initial  condition  should  be 
taken  as 


T*(x*,  r*,  0)  =  0  (2-28) 


However,  for  the  other  cases,  in  order  to  gain  faster 
convergence,  the  initial  temperature  can  be  assumed  to  have 
a  linear  distribution  in  the  axial  direction,  and  thus  have 
the  form 

T*(x*,  r*,  0)  =  T*h  +  YJlc (2-29) 


where  LL  is  the  dimensionless  length  within  which  the  axial 
linear  temperature  distribution  is  assumed  to  hold,  and  xx* 
is  the  dimensionless  distance  measured  from  the  "origin" 
which  is  chosen  only  for  the  purpose  of  this  linear 
temperature  initialization.  Both  the  LL*  and  the  origin 
selected  depend  on  the  model  considered.  For  Model  1  and 
Model  3,  LL*  is  equal  to  L*,  and  the  origin  is  taken  at  the 
left  end  (Model  1) ,  or  at  the  left  intersection  between  the 
central  pipe  and  the  left  extension  pipe  (Model  3) .  However, 
for  Model  2,  equation  (2-29)  is  valid  only  for  the  right 
half  portion  of  the  pipe  as  LL*  is  taken  to  be  L*/2,  and  the 
origin  is  chosen  at  the  middle  section  of  the  pipe.  The 
initial  temperature  of  the  left  half  pipe  can  then  be  found 
from  the  plane  symmetric  condition  about  the  origin  cross- 
section. 
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Calculation  of  Tidal  Displacement 
An  important  quantity  encountered  in  the  study  of  the 
enhanced  heat  transfer  process  in  oscillatory  pipe  flow  is 
the  tidal  displacement,  which  is  usually  required  to  be 
smaller  than  one  half  of  the  total  pipe  length  in  order  to 
avoid  any  convective  mass  exchange  occurring  between  the  two 
reservoirs.  It  has  already  been  defined  in  the  introduction 
chapter  (1-6)  and  in  the  present  computation  the  dimensional 
tidal  displacement  takes  the  form 


AX  =  ^^2   j   x(r,7r/a,)rdr  (2-30) 


where  R^  is  the  inner  radius  of  the  pipe  and  x:(i^fT/w)  is  the 
Lagrangian  displacement  of  the  fluid  elements  located  along 
a  radius  within  the  capillary  tube  at  t  =  tt/w.  It  is 
assumed  these  elements  are  initially  lined  up  at  axial 
position  X  =  L/2  (Model  1  and  Model  2)  or  x  =  5.5L  (Model 
3),  half  way  between  the  tube  ends.  Numerically,  the 
dimensional  Lagrangian  displacement  at  time  t  is  computed 
via  the  equation 


•,t)  =  J   U(r, 


x(r,t)  =     U(r,t)dt  (2-31) 

•'o 

It  is  obvious  that  the  dimensional  Lagrangian  displacement 
x(r,t)  is  a  function  of  both  time  t  and  the  radial  position 
r.   Note  that,  since  the  existence  of  phase  lags  between  the 


H. 


36 
stimulating  pressure  gradient  and  the  displacements  vary  for 
different  Wormersley  numbers,  in  actual  calculations,  the 
tidal  displacement  is  equal  to  the  sum  of  the  absolute 
maximum  cross-section  averaged  Lagrangian  displacement  and 
the   absolute   minimum   cross-section   averaged   Lagrangian 
displacement   within   a   period.     With   the   same   non- 
dimensionalizaton  procedure  used  in  the  previous  section  for 
the  governing  equations,  the  dimensionless  Lagrangian  and 
tidal  displacement  can  be  written  as 


AX*  =  2  J   X*{r*,;r/a,)r*dr*  (2-32) 


The  nondimensional  Lagrangian  displacement  within  the  period 
can  then  be  computed  by 


t* 
X*    (r*,t*)  =  f  U*(r*,t*)dt* 

where 


*      ^ 


(2-33) 


^     (A/c^) 
and 


AX*  = 


U 


(A/w^) 

It  is  pointed  out  that  this  dimensionless  tidal  displacement 
AX*  differs  from  AX  defined  by  equation  (1-7)  by  a  constant 
1/2. 
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Effective  Heat  Flu-x 
The  time  averaged  total  effective  axial  heat  flow 
within  the  pipe  has  the  form 

^total  =  *==^H„   dt    UTrdr  (2-34) 


w 


Where  c  is  the  specific  heat,  p    is  the  density  and  again, 
is  the  oscillatory  frequency.   The  per  unit  area  effective 
heat  flow  (termed  heat  flux)  can  then  be  written  as 

,  ^  ^total 
^      ttRi^  (2-35) 


Equation  (2-34)  follows  from  the  fact  that  pcUT  is  the 
convective  heat  flux.  The  dimensionless  total  effective 
heat  flow  can  be  written  as 

^  r2n  pRl* 

Q  total  =  ^'"^  Jo  '^^*  Jo   U*T*r*dr*  (2-36) 

and  the  dimensionless  heat  flux  by 


Q* 
,*  ^     total 

^  nR*'^  (2-37) 


CHAPTER  III 
NUMERICAL  TECHNIQUES  EMPLOYED 

Equations  (2-5),  (2-6)  and  (2-7)  which  were  derived 
in  Chapter  II  can  not  generally  be  solved  analytically  in 
any  simple  manner.  Therefore,  it  is  necessary  to  seek  a 
numerical  solution  approach  to  the  problem. 

As  is  well  known,  in  order  to  numerically  solve  a  set 
of  simultaneous  governing  differential  equations  more 
accurately  and  efficiently,  an  optimized  grid  system  plays  a 
very  important  role.  For  the  present  purpose  the  grid  net- 
works were  generated  in  the  following  way:  within  the  tube, 
a  non-uniform  grid  (see  Fig.  3-1)  is  clustered  along  the 
radial  direction  in  the  vicinity  of  the  inner  surface  of  the 
pipe  where  the  larger  gradients  of  velocity  and  temperature 
are  expected  to  be  present,  while  along  the  x-axis  the  grid 
is  distributed  uniformly  except  for  those  calculations 
dealing  with  model  3. 

The  solution  process  was  carried  out  in  the 
computational  plane  rather  than  directly  in  the  physical 
plane.  Thus,  a  transformation  which  converts  the  governing 
equations  as  well  as  the  boundary  condition  plays  an 
essential  role. 
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1) .  For  Model  1 
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Fig  3-1  Grid  Systems  Used  in  the  Numerical  Simulations 
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The  second-order  implicit  unconditionally  stable 
Crank-Nicolson  scheme  and  the  ADI  method  [2,  31]  are 
employed  to  break  up  the  transformed  equations  into  finite 
difference  form.  A  computer  code  designated  ETP  (for 
Enhanced  Thermal  Pumping)  has  been  developed  for  obtaining 
the  desired  results.  Details  of  this  procedure  will  be 
described  in  this  chapter. 

Transformation 
The  numerical  solution  of  a  system  of  partial 
differential  equations  can  be  greatly  simplified  by  a  well- 
constructed  grid.  It  is  also  true  that  an  improper  choice 
of  grid  point  locations  can  lead  to  an  apparent  instability 
or  lack  of  convergence.  Early  work  using  finite  difference 
methods  was  restricted  to  some  simple  problems  which  can  be 
numerically  solved  in  the  physical  domain.  As  experience 
was  gained,  general  mappings  were  used  to  transform  the 
physical  plane  into  the  computational  domain  as  well  as  the 
governing  differential  equations  [31].  Such  a  grid 
transformation  technique  brought  to  the  numerical  simulation 
numerous  advantages,  and  the  computational  work  was  no 
longer  restricted  to  a  few  simple  geometries.  A  body-fitted 
non-uniform  grid  in  the  physical  domain  could  be  used,  which 
retains  the  uniform  spaced  grid  system  features  in  the 
computational  domain  [33,  34,  35].  However,  several 
requirements  must  be  placed  on  the  transformation:  first, 
the  mapping  must  be  one  to  one,  second,  the  grid  lines  must 
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be  smooth  and  have  small  skewness  in  the  physical  domain, 
third,  grid  node  point  spacing  should  be  small  where  large 
numerical  errors  are  expected  (i.e.,  large  solution  gradient 
regions)  in  the  physical  domain  (Fig.  3-2) . 

In  present  study,  the  grid  system  in  the  physical 
space  is  numerically  generated  by  solving  an  algebraic 
equation  which  clusters  the  grid  lines  in  the  region  where 
large  gradients  of  the  physical  quantities  are  expected  so 
as  to  gain  higher  resolution  of  these  physical  quantities. 
Because  of  the  simple  pipe  geometry,  the  requirements  of 
grid  smoothness  and  orthogonality  are  not  a  serious  problem. 
For  the  present  purpose,  the  transformation  is  a  simple  one 
which  transforms  a  non-uniform  grid  network  in  the  physical 
plane  onto  a  uniform  grid  system  in  the  computational  plane. 

X*  =  x*(0 

r*  =  r*(,)  ^3_^j 

t*  =  t*(r) 

The  inverse  transformation  can  be  found  as 

?  =  ^  (X*) 

n   =  r,    (r*)  ^3_2j 

r  =  r  (t*) 

Where  x  ,  r*  are  the  dimensionless  coordinates  and  t*  the 
dimensionless  time  in  the  physical  domain,  ^ ,  r,  are  the 
transformed  coordinates  and  r  is  the  transformed  time  in  the 
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computational   domain.     With   this   transformation,   the 
transformed  governing  equations  have  the  form: 


au*  a2u*       au* 


-J7-  =    f  (^)  +  a(,)^^p—  +  b(,)-^^—  (3-3) 

~Tr-  =  pi^r  ^  ^'"^?^  ^  p^^T"  ^  ^'T^  (3-4) 

and 

dT*   ^  dT*  a^T*  dT*  a^T* 

where 

a  t* 

1   at*    1 


a(»?)  = 


2  dr      (ar*/av)''^  (3-7) 


2   ar  [  r-^iar-^/dr,)     ~  (ar-'/ar,)'^  ~d^\ 


and 


-  =  — ( 


at*  ^   CU*  1  g2^* 


at*     1 


P2  = 


P3  = 


P4  = 


(3-8) 


dr     [  ax*/a^   "  2Pr(axV3^)^  "^  ae?  J      (^'^^ 


ar   2Pr(ax*7ay2y  (3-10) 

dr      2Pr  [r'^CarVar,)  "  (arVar?)^   3^2  J        (3-11) 
at*      1 


dr       2Pr  (arVa»7)2  (3-12) 
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at*   1  ,   -1    d^x*  . 


at*   1 


^2  -  dr      1¥^  JdxVdTP  (3-14) 


at*   1  ,    1         1    a2i-* 


'r^^  [r*(ar 


^3  -  ar  2Prw  [r'^carVar?)  ■  (ar*/a„)^  ?^J       (^-is) 

at*   1     1 
^^  =  ~J7~  I^  (ar*/a,)^  (3-16) 

As  mentioned  in  chapter  II,  equations  (3-3),  (3-4)  and  (3-5) 
are  a  set  of  second-order  parabolic  partial  difference 
equations  in  cylindrical  coordinates.  In  addition,  since 
the  oscillating  flow  is  considered  as  incompressible  in  the 
present  studies,  the  momentum  equation  (3-3)  can  then  be 
independently  solved  at  each  time  step.  As  a  result,  the 
time-dependent  update  velocity  found  can  then  be  substituted 
into  the  heat  conduction  equation  (3-4)  as  a  known  quantity 
at  the  same  time  step  level.  Eventually,  the  coupled  heat 
conduction  equations  (3-4)  and  (3-5)  are  solved 
simultaneously  to  obtain  the  temperature  distribution  both 
in  the  fluid  and  in  the  wall  at  any  time.  To  best  solve 
this  set  of  equations  in  terms  of  accuracy  and  efficiency, 
the  proper  choice  of  numerical  technique  and  grid  net-work 
is  dictated  by  an  understanding  of  the  physical  aspects  of 
the  problem. 
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The  same  transformation  should  be  also  applied  on  all 
the  boundary  conditions  proposed  in  the  three  different 
models.  For  the  temperature  compatibility  conditions  along 
the  interface  between  the  fluid  and  the  wall  (2-25)  one  has 
the  following  transformed  forms: 


aT"    1 

dri        (ar*/dr,)     ^f 


fluid 


dT  1 

dn      {dr*/dr,)    ^^ 


wall 


(3-17) 


or 


dT 


=     K; 


dT" 


^^    'flow   '^^  ^^  'wall 


(3-18) 


where 


Ka  = 


kv  dr*/dt,^ 


pipe 


kf  dr*/dri  I   TT 
•^    '  '  I  wall 


(3-19) 


To  make  the  subsequent  form  of  the  corresponding 
finite  difference  governing  equations  less  cumbersome,  the 
superscript  "  *  "  will  be  dropped  from  the  variables  T*,  r*, 
t  ,  and  X  ,  and  in  addition,  U  will  be  replaced  by  V. 

In  the  process  of  deriving  the  finite  difference 
governing  system  equations,  the  second-order  central 
differences  in  the  domain  and  forward  or  backward 
differences  along  the  boundaries  or  interface  of  the  fluid 
and  wall  have  been  employed  at  each  nodal  point. 
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Crank-Nicolson  Method  for  Momentum  Ecmatinn 
The  second-order  accurate  Crank-Nicolson  method  is 
quite  well  known  and  widely  used.  It  is  an  unconditional 
stable,  implicit  scheme  for  solving  the  parabolic  types  of 
partial  difference  equations.  when  the  Crank-Nicolson 
method  is  applied  to  equation  (3-3),  the  finite  difference 
algorithm  at  a  typical  node  k  in  the  radial  direction  and  at 
the  time  step  n  assumes  the  simpler  form: 

k  =  1/  2,  km 


(3-20) 


where 


r 


=   2  (  *"^'*   ^"  ]  -  (  Bj^V,^.^  +  E^V^-   A^V^^^  ]"   (3-21, 


1 

^k  =  -  — (  2  a   +  b  )3^  (3_22j 


B,^  =  — (2  a  +  b  )^  ^3_23, 


°k  =  1  -^  -k  (3-24) 


''k  -  ^  -  ^k  (3-25) 
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The  boundary  conditions  along  the  interface  (i.e., 
the  inner  surface  of  the  wall,  k  =  kmid)  and  the  axis  (k  = 
1)  then  become 


^  |k  =  kmid"  ° 

(3-26) 


av 


dr, 


k  =  1=  0 


The  initial  conditions  of  the  velocity  at  all  nodal 
points  is  taken  to  be  zero. 


\  n  =  1  =  0 


k-1,  2,  3,  kmid 

(3-27) 

Equation  (3-20)  associated  with  the  boundary  conditions  (3- 
26)  can  then  be  written  in  the  matrix  form.  This  yields  a 
set  of  linear  system  algebraic  equations  which  can  be  solved 
in  terms  of  the  nodal  values  of  velocity  in  the  capillary 
tube  by  using  either  an  iterative  method  or  a  Thomas 
algorithm  at  each  time  step.   The  explicit  form  is 
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^kd-l^d-Ad-l 


B 


kd 


kd 


n 

^1 

n+1 

^1 

^2 

^2 

^3 

= 

S 

V 
^kd-1 

\d-l 

J 

.  ""kdj 

.  ^kd  . 

n 


where 


(3-28) 


A,  = 


H   -^  ^1 


kd  =  kmid-1 

API  Method  for  Axisymmetric  Heat  Ecmations 
The  governing  PDEs  (3-3),  (3-4)  and  (3-5)  are  all  of 
the  second-order  parabolic  type.  Thus  it  might  be  suggested 
that  the  Crank-Nicolson  scheme  used  in  solving  the  momentum 
equation  (3-3)  can  also  be  applied  to  the  axisymmetric  heat 
equations  (3-4)  and  (3-5),  and  one  can  then  take  advantage 
of  the  tridiagonal  matrix  form  while  using  this 
unconditionally  stable  technique.  However,  when  attempting 
to  use  such  a  formulation,  one  immediately  finds  that  the 
resulting  system  of  linear  algebraic  equations  is  no  longer 
of  the  tridiagonal  type  (3-20)  but  rather  a  non-tridiagonal 
matrix  system  requiring  substantial  CPU  time  to  solve.  This 
difficulty  can  be  avoided  by  applying  the  unconditionally 
stable   Alternating   Direction   Implicit   method   (ADI) , 
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developed  by  Peaceman  and  Rachford  and  Douglas  in  1955. 
According  to  this  scheme,  the  entire  solution  process  at 
each  time  step  is  "split"  into  two  portions,  i.e.,  the  first 
half  of  solution  processes  for  k''^^  column  (radial 
direction)  ,  while  the  other  half  processes  for  the  j*^'^  row 
(axial  direction) . 

With  the  ADI  scheme,  second-order  central  differences 
are  used  to  approximate  the  values  of  derivatives  at  each 
nodal  point  in  equations  (3-4),  (3-5).  The  finite 
difference  algorithm  for  those  equations  during  the  first 
half  of  each  time  step  for  the  j"^^  column  are  then 

I  BP .  .  1  T  .  .  ,  ,  +  DP .  .  ,  T  .  .  ,  +  AP  .  .   T  .  .      1  "■•"  V2  ^  p„n 
I   ::,k  DD,k-l   ''^DD,k^jj,k  ^   :j,k^jj,k+l  J  jj,k 

k  =  1,  2,  ... .kmid 
(3-29) 


and 


I  BW  .  .  ,  T  .  .  ,  ,  +  DW  .  .  ,  T  •  .  ,  +  AW  .  .   T  .  .      1  ""^  V2  ^  „„n 


k  =  kmid+1, . . .kmax 
(3-30) 


where  the  subscript  jj  is  used  to  emphasize  the  specific 
column  currently  to  be  computed.  It  can  be  seen  that  the 
set  of  difference  equations  now  is  in  the  tridiagonal  form 
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since  the  right-hand-side  terms  in  the  equations  (3-29) ,  (3- 
30)  contain  only  known  values  from  the  previous  results  and 
the  boundary  conditions.   These  values  can  be  computed  by 

^^jj,k  =  (-^^jj,kTjj-l,k  ^   ^^jj,kTjj,k  -  ^^jj,k^jj.l,k]" 


(3-31) 


and 


wx 


n 


jj/k 


=  ( 


-CW .  .  ,  T .  .  ,  ,  +  EW .  .  ,  T .  .  , 


-  FW .  .  ,  T  .  .  ,  ,  ,  1  " 


(3-32) 

The  computational  algorithm  is  implemented  column  by 
column  and  the  unknown  value  (Tj,]^)"  can  then  be  solved  by 
either  an  iterative  or  a  direct  method.  In  order  to  do 
this,  equation  (3-29)  needs  to  be  assembled  into  the 
following  matrix  form. 


)*#^ 


i  i 


PD^   PA^ 

PB^   PD^   PA^ 


PB3   PD3   PA3 


P\d-lPDkd-lP^kd-l 
^^kd   P°kd 


n 

'   ^1 

n+1/2 

'   P^l 

^2 

PX^ 

^3 

= 

PX3 

'^kd-1 

P^kd-1 

J 

.   ^M. 

J 

n 


j  =  1,  2,  3 


3  max 
(3-33) 


51 
Similarly,  equation  (3-30)  can  be  assembled  into  the  matrix 
form  as  seen  in  (3-34) . 


^°kd+l  W^kd+1 
W^kd+2  ^°kd+2  W^kd+2 


WB,     WD,     WA, 
km     km     km 

WB,     WD, 
kmax   kmax 


n 

T 

kd+1 

n+1/2 

^^kd+1 

'^kd+2 

^^kd+2 

T 

■^km 

^^km-1 

J 

T 
kmax 

J 

^^kmax 

n 


j  =  1,  2,  3 


jmax 


(3-34) 


tf  ■■  ■ 


It  should  be  pointed  out  that  the  coefficients  and/or 
the  right  hand  side  terms  in  the  row  marked  with  symbol  "*" 
in  those  matrixes  needs  to  be  properly  modified  as  well  as 
the  limit  of  either  subscript  j  or  subscript  k  according  to 
the  different  boundary  conditions  in  the  various  models 
considered.  For  example,  along  the  axis,  the  symmetry 
condition  requires  that  ST/ar  =  0  and  this  could  be 
accomplished  numerically  by  equating  the  temperature  Tq  and 
T2  and  by  employing  a  new  combining  coefficient  of  PA^*  = 
PAi  +  PBi  rather  than  the  original  PAi  in  the  first  matrix 
for  evaluating  the  temperature  distribution  within  the  pipe. 
In  addition,  if  the  temperature  distributions  along  certain 
parts  of  the  outer  surface  of  the  wall  are  given,  then  the 
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limit  of  subscript  k  will  be  ended  with  k  =  kmax-1  rather 
than  k  =  kmax  shown  in  equation  (3-34) .   The  same  arguments 
are  also  applicable  for  another  subscript  j  in  the  axial 
direction. 

By  using  a  similar  procedures  for  the  second  half  of 
each  time  step,  the  finite  difference  algorithm  for  the  k'^'^ 
row  then  becomes 

j  =  1,  2,  jmax 

(3-35) 

(^«j,kkTj-l,kk  ^  ^^j,kk^j,kk  ^  ^^j,kkTj+l,kk)"^'  =  W^kk! 

j  =  1,  2,  jmax 

(3-36) 

where  the  right  hand  terms  can  be  computed  by 

^^j,kk  =  ["^^j,kk'^j,kk-l  +  "Pj,kk'^j,kk'  ^^j,kk'^j,kk+lj''^^'^^ 

(3-37) 

^<ik'=  [-c«j,kkTj,kk-i  +  «^j,kk^j,kk-  ^Wj,kkTj,kk+i]"^'''' 

(3-38) 


B-v., 
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Similarly,  to  emphasize  that  the  current  calculation  is  in 
the  ]iP^    row,  we  use  the  symbol  kk  rather  than  k  in  above 
formula  and  the  matrix  form  of  equation  (3-35)  for  the  k^^ 
row  can  then  be  written  as 


PG2   PF2 

PC3   PG3   PF3 


PC.   PG.   PF^ 
4     4     4 


P^jm-1  P^jm-1  P^jm-1 


PC. 
:m 


PG. 
Dm 


n+1/2 

'  ^2 

P^2 

^3 

PY3 

^4 

PY4 

T. 
3m-l 

P^jm-1 

k 

T. 

k 

Dm  J 

n+1/2 


k  =  2,  3 


,  kmid 


(3-39) 


and  the  matrix  form  of  equation  (3-36)  becomes 


WG^   WF^ 

WC3   WG3   WF3 


WC.   WG.   WF^ 
4     4     4 


WC .     WG .     WF . 

WC .     DG . 
jm     jra 


n+1/2 

^2 

N+1 

WY2 

^3 

WY3 

^4 

WY. 
4 

• 

T. 

Dm-1 

k 

T. 
L   Dm  ^ 

k 

Dm  J 

'  N+1/2 


k  =  kmid+1,  kmid+2 km 


(3-40) 
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As  indicated  above,  the  symbol  "  *  "  shows  that  the 
coefficients  in  that  row  as  well  as  the  right-hand-side 
terms  need  to  be  properly  modified  corresponding  to  the 
different  boundary  conditions.  The  matrix  terms  in 
equations  (3-29),  (3-30),  (3-31)  and  (3-32),  (3-35),  (3-36), 
(3-37)  and  (3-38)  can  be  computed  by 


^^j,k  =  -(  ^P3  +  P4) 


j,k 


^Pj,k=  (^  -  2*P2).^j^ 


^^j,k  =  (— P3  -  P4) 


j,k 


^Pj,k=  -("i-Pl  +  P2) 


j,k 


^Pj,k  =  (1-  PI  -  P2) 


j/k 


^^j,k=  (1  +  2*P2).^,^ 


DP 


j,k=  (1  ^  2*^^)j,k 


«Pj,k=  (1  -  2*P4)j,k 


and 


(3-41) 


AW.    =  -(-^W3  +  P4)  . 


j,k 


j,k 


^^j,k=  (1  -  2*P2).^^ 


BW     =  (-7-W3  -  W4) 


J/k 


^^j,k  =  -(^-Wl  +  W2).^. 


^^j,k  =  (^-  Wl  -  W2) 


j,k 


^^j,k=  (1  +  2*W2).^j^ 


DW. 


j,k=  (1  +  2*W4).^3^ 


«^j,k=  (^  -2*W4)j^,, 


€^1-  ■^•■" 


'^'■t-i' "'"^^  - 


(3-42) 


^ 


(3-43) 
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In  order  to  rewrite  the  equation  in  a  simple  form  we  define 
the  following  linear  operator  : 

LX  [P]j^k  =  t^^'  ^^'  ^^^j,k 

LX  [W]j^j^  =  [CW,  GW,  FW]j^j^ 
and  column  vectors 

^'^^>j,k=  <Tj,k-l'  Tj,k'  Tj,k.l  > 
t^^>j,k=  ^Tj-l,k'  ^j,k'  Tj+l,k  > 


(3-44) 


Equations  (3-29)  which  are  in  the  radial  direction  can  be 
simply  written  as 


n       n+1/2     n 
LY  [P]     {TY}       =  PX  (3-45) 


j/k     j,k       j,k 


n       n+l/2     n 
LY  [W]     (TY)       =  WX  (3-46) 

j,k     j,k       j,k 


J  =  1,  2, 
k  =  1,  2, 
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Equations   (3-33)   and   (3-34),   which   are   in   the   axial 
direction,  then  assume  the  following  simple  forms: 


n+1/2     n+1     n+1/2 
LX  [P]       {TX}     =  PY  ^3.7. 

3,k        j,k      j,k  ^^  ^^^ 


n+1/2     n+1      n+1/2 
LX  [W]       (TX)      =  WY  ,3.4g. 

j'k        j,k       j,k  ^^  ^^^ 


:  =  1,  2, 
k  =  1,  2, 

where  the  right  side  terms  can  be  estimated 


PX_    =  (-CP,  EP,  -FP)    {TX} 
^'^  j,k     j,k 


WX^    =  (-CW,  EW,  -FW)    (TX) 

^'^  j,k    j,k 


PY_    =  (-BP,  HP,  -AP)    (TY) 
"^'^  j,k     j,k 


WY     =  (-BW,  HW,  -AW)    (TY) 
^'^  j,k     j,k 


by 


(3-49) 
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It  is  seen,  as  the  result  of  the  ADI  "splitting" 
procedure  which  has  been  employed  in  the  algorithm 
associated  with  different  boundary  conditions  for  the 
various  models,  that  only  a  tridiagonal  system  of  linear 
algebraic  equations  needs  to  be  solved  (i.e.,  during  step  1, 
the  coupled  tridiagonal  matrix  (3-33),  (3-34)  are  solved  for 
each  jth  column  of  the  grid  points,  while  during  step  2,  the 
coupled  tridiagonal  matrix  (3-39)  and  (3-40)  are  then  solved 
for  each  k^h  j.^^  ^^  g^^^^^  points)  . 

Once  the  periodic  steady  solution  has  been  obtained, 
we  can  calculate  both  the  tidal  displacement  and  the  heat 
flux  at  different  locations  within  the  pipe.   The  numerical 
technique  used  here  to  integrate  equations  (2-38),  (2-39) 
and  (2-41)  for  evaluating  the  tidal  displacement  and  the 
heat  flux  at  each  specified  location  can  be  obtained  either 
by  the  Trapezoidal  rule  (with  end  correction)  or  Simpson's 
rule.    Both  numerical  approaches  are  essentially  fourth- 
order  methods. 

Convergence  CrJ-hAT-ia 
As  is  known,  once  the  calculation  work  has  started, 
the  time-matching  process  will  be  in  the  loop  forever  unless 
a  criterion  can  be  derived  that  indicates  when  the  goal  of 
the  current  computing  work  has  been  reached  and  further 
solution-matching  processes  do  not  produce  significant 
increases  in  accuracy.  Such  a  criterion  depends  on  the  : 
purpose  of  the  calculation.   If  one's  goal  is  to  study  the 
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transient  process,  one  can  indicate  a  time  limit  to  pause  or 
quit  the  current  computing  work. 

However,  if  one  seeks  a  final  periodic  state 
solution,  one  has  to  develop  a  criterion  to  test  if  the 
solution  can  be  considered  acceptable. 

It  seems  that  the  temperature  is  a  good  measure  of 
the  accuracy  of  the  overall  solution  process  so  that  the 
most  efficient  way  is  to  apply  the  convergence  test  on  the 
temperature  rather  than  on  the  velocity.   For  the  present 
purpose,  two  convergence  criteria  were  alternatively  used. 
The  choice  of  the  criterion  depends  on  what  was  the  main 
goal  in  the  particular  study  case,   if  the  main  effort  is  to 
observe  the  distribution  of  the  quasi-steady  temperature  at 
any  phase  angle  .t  within  a  period,  the  testing  is  done  by 
comparing  the  temperature  residual,  i.e.,  by  inspecting  the 
averaged  temperature  difference  of  each  nodal  point  at  thee 
same  u,t  between  adjacent  periods.   This  can  be  written  as 


Resl  = 


I  S  (Til  -  T22  )?    ]V2 


(JMAX)  (KMAX)  ^  ^■'■ 


3  =  1/  2,  JMAX 

k  =  1/  2,  KMAX 


(3-50) 
where  Tll.^^^  is  the  value  of  the  temperature  of  node  (j,k) 
at  a,t  in  the  current  period,  T22.^j^  is  the  value  of  the 


..'■>■' 
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temperature  of  node  (j,k)  at  the  same  .t  in  the  previous 
period,  and  e^  is  the  convergence  parameter. 

If  the  goal  of  the  investigation  is  to  examine  the 
effective  heat  transfer,  the  convergence  criterion  is 
established  by  computing  the  residual 


f  S  (  ^1  -  <i>2    )^  ]^/^ 


Res2  =  -  ,      ^  ^ 

JSEC  '  "2  ^'-^^^ 


J  =  1/2, JSEC 

where  JSEC  is  the  number  of  cross-sections  where  the  axial 
heat  flux  was  examined,  ^2  is  the  convergence  parameter,  and 
^1,  ^2  are  the  cross-section  averaged  heat  flux  in  the 
current  period  and  previous  period,  respectively.  The 
summation  is  carried  out  over  all  sections  (JSEC)  where  the 
heat  flux  is  computed. 

The  value  of  the  convergence  parameter  can  be 
determined  by  a  balance  of  the  acceptable  solution  accuracy 
and  the  cost  of  CPU  time.  m  the  current  study  both  el  and 
e2  were  taken  between  0.001-0.01  as  the  residual  of  the 
temperature  or  the  heat  fluxes  for  acceptable  value. 

Grid  Generation 

As  is  well  known,  the  solution  accuracy  and  efficiency 
in  a  large  degree  depends  on  the  grid  system  used  in  the 
numerical  calculation.  a  good  grid  is  characterized  by 
small  skewness,   high  smoothness  and  capability  of  high 
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resolution  in  the  large  gradient  regions  in  the  physical 
plane.  It  has  been  shown  that  rapid  changes  of  grid  size 
and  highly  skewed  grids  can  result  in  undesirable  errors 
[25].  The  success  of  a  numerical  simulation  of  a  complex 
thermal  fluid  dynamics  problem  does  depend  strongly  upon  the 
grid  system  used  in  the  computation. 

In  the  present  study,  the  grid  mesh  was  generated  with 
the  emphasis  on  the  high  resolution  capability  of  its  simple 
geometric  boundaries.  The  technique  used  here  is  one  of  the 
algebraic  schemes  which  cluster  the  grid  lines  near  the 
region  desired  [34],  namely. 


AS,  =  AS  (1  +  e)  ^'-^  (3-52) 

Jc     o 


where  AS  is  the  minimum  specified  grid  spacing  next  to  the 
wall  or  to  some  inner  interfaces  within  the  tube.  The 
parameter  e  is  determined  by  a  Newton-Raphson  iteration 
process  so  that  the  sum  of  the  above  increments  matches  the 
known  arc  length  between  k  =  1  and  k  =  kmax.  The  grid 
networks  used  in  Model  1  and  Model  2  are  uniformly 
distributed  along  the  axial  direction  in  both  the  fluid  and 
the  wall;  however,  in  Model  3,  most  grid  points  are 
clustered  near  the  central  portion  of  the  pipe  so  as  to  gain 
higher  resolution  in  the  region  of  interest,  while  along  the 
radial  direction  a  non-uniform  grid  was  employed  in  all 
three  models.   The  minimum  specified  grid  spacing  aSq  used 
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depends  on  the  boundary  layer  thickness  S,  namely,  the 
kinematic  viscosity  of  working  fluid  and  the  oscillating 
frequency.  In  the  present  study,  it  was  found  that  a  good 
choice  of  this  value  is  aSq  =  0.055  for  laminar  flow  cases 
with  a  total  of  about  15-20  nodal  points  distributed  along 
the  radius. 


V  ■  ■(■ 


CHAPTER  IV 
NUNERICAL  RESULTS  AND  DISCUSSION 


The  problem  of  time-dependent  enhanced  heat 
conduction  subjected  to  sinusoidal  oscillations  can  now  be 
solved  numerically  for  the  boundary  conditions  appropriate 
to  a  long  capillary  tube  according  to  the  various  models 
described  in  previous  chapters.  The  computational  tasks 
fall  into  two  catagories:  the  first  part  is  a  numerical 
study  concerned  with  the  characteristics  of  the  oscillatory 
pipe  flow,  which  includes  an  investigation  of  the  velocity 
profiles,  the  Lagrangian  and  the  tidal  displacement 
trajectories.  The  second  part  of  this  study  includes  a 
thorough  investigation  of  the  final  periodic  state 
temperature  build-up  process  in  oscillating  pipe  flow,  the 
periodic  temperature  distribution  in  the  pipe  and  wall,  and 
a  study  of  the  relationship  between  the  enhanced  axial  heat 
flux  and  the  tidal  displacement.  It  also  includes  an 
investigation  of  the  tuning  effect  and  a  comparison  of  the 
enhanced  axial  heat  transfer  with  the  corresponding  pure 
molecular  axial  heat  conduction  as  well  as  the  investigation 
of  the  influences  of  the  pipe  radius  and  the  wall  thickness 
on  the  enhanced  axial  heat  transfer. 


62 


63 
Properties  of  oscillating  laminar  pipe  flow  have  been 
analytically  discussed  by  S.  Uchida  [36]  for  several 
different  Wormersley  numbers,  and  these  results  offer  a 
useful  reference  for  comparing  with  the  present  numerical 
studies.  However,  in  the  area  of  thermal  fields  associated 
with  enhanced  thermal  pumping  there  is  little  information 
available  for  comparison,  except  for  some  recent  results  of 
Kaviany  [13]. 

The  present  computational  work  was  carried  out  on  the 
Vax-11/750   computer   in   the   Department   of   Aerospace 
Engineering,  Mechanics,  and  Engineering  Science,  and  on  the 
Vax-11/780  in  the  Center  for  Instructional  and  Research 
Computing  Activities,  University  of  Florida,  Gainesville, 
Florida.   Note  that  it  is  very  time-consuming  to  build-up  a 
final  periodic  state,   for  example,   if  the  grid  size  is 
101x22  and  the  time  steps  are  2000  per  period,  it  takes  7.5 
minutes  (CPU)  with  the  VAX-11/780  machine  or  almost  1  hour 
(CPU)  with  the  VAX-11/750  machine  to  run  only  one  period. 
It  usually  takes  20-30  periods  to  reach  the  final  periodic 
state  solution. 

The  process  of  selecting  model  size  (i.e.,  total 
nodal  points)  is  a  synthetic  balance  among  the  storage 
requirement,  the  solution  accuracy,  and  the  cost  of  CPU  time 
for  finding  an  acceptable  solution.  Such  a  selection  allows 
each  case  to  be  solved  with  a  minimum  of  expense  in 
computing  operation  thereby  making  it  possible  to  do  the 
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large  number  of  runs  needed  to  obtain  enough  data  points  to 
plot  curves  of  the  desired  parameters. 

The  results  discussed  in  this  chapter  are  presented 
for   the   purpose   of   illustrating   the   effects   to   be 
encountered  in  working  with  the  enhanced  thermal  pump. 
These  results  should  be  effective  in  gaining  insight  into 
some  interesting  features  of  this  enhanced  heat  transfer 
process.   Since  the  analytic  approximations  [15,  16,  18,  22] 
are  effective  in  describing  the  flow  and  the  heat  transfer 
aspects  and  to  give  considerable  insight  into  the  problem, 
they  are  often  used  to  compare  with  the  numerical  results 
and  should  give  good  agreement  when  applicable. 
Part  1.   Oscillatory  Pipe  Flow  Features 
In  order  to  better  understand  the  mechanism  of  the 
enhanced  thermal  pump,   it  is  necessary  to  examine  the 
mechanical   features  of  the  oscillating  pipe  flow.    As 
mentioned  in  previous  chapters,  the  present  interest  in  the 
enhanced  thermal  pumping  is  confined  to  an  investigation  of 
the  central  part  of  the  slender  pipe,  thus  the  flow  field 
can  be  well  approached  by  a  1-D  time-dependent  laminar  model 
which  neglects  the  ends  effects. 
Velocitv  Profiles 

Fig.   4-1   shows   the   numerically   computed   time- 
dependent  velocity  profiles  at  different  phase  angles  of  the 
exciting  pressure  when  Wormersley  numbers  are  equal  to  1, 
^i-      10,  100,  and  1000,  respectively.   It  can  be  clearly  seen 
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that  the  velocity  profile  at  a  =  1  presents  a  quasi- 
parabolic  shape  at  any  instant  within  a  period  and  is  in 
phase  with  the  stimulating  pressure  gradient.  However,  at 
higher  frequency,  for  example,  a  =  100  and  1000,  the 
velocity  profiles  can  be  clearly  separated  into  two  regions: 
in  the  vicinity  of  the  wall  the  flow  shows  a  typical  thin 
boundary  layer,  while  in  the  region  far  away  from  the  wall 
the  fluid  moves  as  if  it  were  frictionless  slug  flow.  In 
fact,  within  this  core  region  the  velocity  distribution  is 
independent  of  the  distance  from  the  wall  [29].  It  can  be 
also  seen  that  the  phases  of  the  velocity  profiles  at  higher 
frequencies  cases  are  shifted  about  jr/2  with  respect  to  the 
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Fig  4-2  Velocity  Prof ile (a=5  )  [29] 
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Fig   4-3      Magnified  View   of  Velocity   Profile   Near  Wall 
(Wormersley  Number  q=10,    H2O,    5=0.014cin) 
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stimulating   pressure   gradient.     At   the   intermediate 
frequency  case  (a  =  10  of  Fig.  4-1) ,  the  slug  flow  boundary 
is  not  so  evident,  but  one  can  still  note  a  boundary  layer 
near  the  wall.   The  same  pattern  of  the  velocity  profiles 
associated  with  moderate  Wormersley  numbers  can  also  be 
found  in  the  reference  [36]  where  the  velocity  profile  at 
Wormersley  number  a  less  than  10  was  presented  (Fig.  4-2). 
In  order  to  better  see  the  variation  of  the  time-dependent 
velocity  profile  within  the  boundary  layer,  a  set  of  closer 
view  of  the  velocity  profile  with  respect  to  Wormersley 
number  at  a  =  10  during  phase  intervals  wt  =  30°  is  plotted 
in  Fig.  4-3. 

It  should  be  emphasized  that  the  solutions  shown  here 
are  under  the  assumption  that  the  secondary  velocity  in  the 
radial  direction  is  negligible  compared  to  the  axial 
velocity  component,  namely,  the  non-linear  inertial  terms 
are  not  considered  in  the  governing  equation  (2-1) .  This  is 
a  reasonable  approximation  for  moderate  oscillatory 
frequency  w  except  near  the  ends  of  the  pipe.  However,  at 
high  Wormersley  numbers,  care  should  also  be  taken  to  avoid 
violating  another  assumption,  namely,  that  of 
incompressibility  and  the  concomitant  condition  that  the 
oscillating  phase  does  not  change  between  the  tube  ends. 
High  Wormersley  number  can  be  obtained  either  by  increasing 
the  oscillating  frequency  w  or  by  increasing  the  pipe 
diameter  for  a  given  fluid.   An  oscillating  flow  can  be 
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considered  incompressible  if  Axw/2  <  0.05C,  where  C  is  the 
speed  of  sound  in  the  fluid.  For  mercury,  C  =  1360  m/sec, 
one  requires  a  =  27.1,  when  Ax  =  20  cm  and  R^  =  0.1  cm.  To 
avoid  a  appreciable  phase  difference  between  the  tube  ends, 
one  requires  that  L/C  «  2k /w.  Both  the  restrictions  are 
met  in  the  examples  to  be  considered  below. 
The  Laqrancfian  Displacements 

An  alternative  interesting  representation  to  the 
oscillatory  velocity  field  are  the  Lagrangian  displacements 
of  the  fluid  elements  at  different  radii  within  the  pipe. 
They  have  been  plotted  in  Figs.  4-4  and  4-5  at  time 
intervals  of  wt  =  30°  for  Wormersley  number  a  =  0.1,  1.0  and 
10.  It  is  noted  that  since  both  pipe  diameter  and  the 
working  fluid  were  fixed  in  this  test,  Figs.  4-4  and  4-5 
represent  the  relationship  between  the  Lagrangian 
displacement  and  the  oscillating  frequency.  The 
trajectories  plotted  in  Figs.  4-4  and  4-5  have  been 
normalized  by  A/w^,  where  A  =  l/p|dp/dx|  is  the  amplitude  of 
the  sinusoidal  pressure  gradient  as  defined  in  equation  (2- 
1)  and  w  is  the  angular  velocity.  Suffice  it  here  to  point 
out  that  for  the  lower  frequency  case  (for  example  a  =  0.1 
and  1.0)  the  Lagrangian  displacement  trajectory  shows  a 
foreseeable  parabolic  pattern  at  any  moment.  Nevertheless, 
the  essential  distinction  between  low  frequency  oscillatory 
pipe  flow  and  steady  Hagen-Poiseuille  flow  is  that  in  the 
former,  the  Lagrangian  displacement  trajectories  as  well  as 
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Fig  4-4    Lagrangian  Displacement  for  a  =  0.1  and  1.0 
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Fig  4-5    Lagrangian  Displacement  at  a  =  10 
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the  velocity  profiles  are  periodic  so  that  the   fluid 
particles  do  not  translate  axially  upon  time  averaging, 
while  in  the  later  case  they  will. 

For  intermediate  Wormersley  number  (a  =  lo) ,  the 
trajectories  of  the  Lagrangian  displacement  departs 
considerably  from  the  standard  parabolic  shape.  This 
phenomenon  can  be  even  more  clearly  seen  in  the  a  =  loo 
case.  Evidently,  the  higher  the  oscillating  frequency,  the 
thinner  the  boundary  layer  (5  =  Jlv/w)  . 
Tidal  Displacpm<:>ni-g 

Figs.  4-6  and  4-7  demonstrate  time  variation  of  the 
cross-section  averaged  dimensionless  Lagrangian  displacement 
at  Wormersley  numbers  in  the  range  from  a  =  0.1  to  50.  The 
tidal  displacement  can  be  obtained  by  summing  the  absolute 
maxima  and  the  absolute  minima  of  these  curves.  The 
corresponding  non-dimensional  tidal  displacements  with 
respect  to  Wormersley  number  from  a  =  0. 1  to  a  =  12  are 
listed  in  table  4-1  below 


Table  4-1   Dimensionless  Tidal  Displacement 
at  Different  Wormersley  Numbers 


a 

0.1 

1.0 

2.0 

3.0 

4.0 

5.0 

AX 

0.00246 

0.14295 

0.67303 

1.20213 

1.41160 

1.49041 

a 

6.0 

7.0 

8.0 

10.0 

12.0 

AX 

1.56294 

1.61879 

1.66046 

1.71900 

1.76698 
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Fig  4-6  Dimensionless   Cross-section  Averaged 

Displacenent  Versus   Time 
a    =    0.1    -    1.0 
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Fig  4-7    Dimensionless  Cross-section  Averaged 

Displacement  Versus  Time 
Q  =  2  -  50 
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It  is  noticed  that  for  very  small  Wormersley  number 
cases  (a  <  0.5),  the  cross-section  averaged  displacement 
varies  like  a  sinusoidal  function  with  respect  to  time  (the 
ordinate) ,  and  as  Wormersley  number  increases,  the  cross- 
section  averaged  displacements  are  no  longer  symmetric  about 
the  ordinate  but  rather  favor  positive  values  of  x«   If  the 
Wormersley  number  is  further  increased,   eventually,   the 
cross-section  averaged  displacement  almost  entirely  lies  on 
the  right  half  plane.   A  similar  feature  can  also  be  seen  in 
Figs.  4-4  and  4-5.    In  fact,  this  just  again  shows  the 
existence  of  the  phase  shift  in  the  oscillating  pipe  flow. 
When  the  Wormersley  number   is   small   the  cross-section 
averaged   displacements   as   well   as   the   Lagrangian 
displacements  show  a  sinusoidal  variation  with  respect  to 
time  and  this  has  7r/2  phase  lag  with  respect  to  the 
stimulating  pressure  gradient  which  is  assumed  to  be  a 
cosine   function  of  wt  with   zero   initial   phase  angle. 
However,   for  higher  Wormersley  number,   the  phase   lags 
increase  to  almost  tt  with  respect  to  the  phase  of  the 
exciting  pressure  gradient.     Note  that  the   Lagrangian 
displacements  in  Figs.  4-4  and  4-5  were  computed  by  lining 
up  all  the  fluid  particles  on  the  plane  x  =  0  at  wt  =  0, 
however,  by  using  non-zero  velocity  at  wt  =  0°  as  shown  in 
Fig.  4-1.   This  implies  that  if  the  phase  of  the  exciting 
pressure  gradient  is  taken  as  the  base  of  the  measurement, 
it  is  generally  not  possible  to  assure  the  same  phase  to  the 
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velocity  and  Lagrangian  displacement  as  well  as  the  tidal 
displacement.  For  example,  at  phase  of  the  stimulating 
pressure  gradient  wt  =  0°  ,  the  corresponding  velocity  phase 
may  be  7r/2  and  the  Lagrangian  displacement  phase  may  be  tt  . 
In  fact,  this  is  just  as  true  for  very  high  Wormersley 
numbers;  however,  the  phase  difference  with  respect  to  the 
exciting  pressure  gradient  phase  is  less  than  the  value 
shown  above  for  low  Wormersley  number.  One  can  well  see 
that  if  the  phase  lags  of  the  Lagrangian  displacement  or  the 
tidal  displacement  is  given,  one  can  certainly  re-draw  the 
diagram  shown  in  Figs.  4-6,  4-7  by  lining  up  the  phase  with 
itself,  and  then  an  exact  symmetric  pattern  of  the  tidal 
displacement  curve  similar  to  that  in  the  very  small 
Wormersley  number  case  can  be  obtained.  Unfortunately,  the 
phase  lags  are  a  function  of  the  Wormersley  number  and  they 
are  not  known  in  advance. 

In  order  to  verify  the  numerical  method  used  in  this 
study  and  to  check  the  ETP  code  developed,  a  comparison  of 
the  computed  dimensionless  tidal  displacement  to  the  one 
using  the  analytic  equations  (1-7)  and  (1-8)  given  by 
Kurzweg  [18]  for  Wormersley  number  varying  from  0.1  to  100 
has  been  calculated  and  plotted  in  Fig.  4-8.  The  solid  line 
shows  the  analytic  solution  obtained  by  use  of  equations  (1- 
7)  and  (1-8) ,  while  the  dashed  line  shows  the  results  with 
the  ETP  code  developed  in  this  study.  The  agreement  is 
quite  good,  particularly  when  the  Wormersley  is  less  than 
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3.0.  However,  at  high  Wormersley  number  the  numerical 
solution  shows  a  very  slight  deviation  from  the  analytic 
solution.  This  deviation  is  believed  due  to  an  inaccurate 
numerical  integration  over  the  cross-section  using 
relatively  large  time  steps  (our  time  steps  per  period  in 
the  calculation  were  between  1000-2500) .  A  comparison  with 
using  10'*  time  steps  per  period  for  Wormersley  number  a  =  10 
was  studied  and  shows  some  improvement.  However,  using  such 
small  time  steps  in  the  present  investigation  is  beyond  the 
capacity  of  the  current  VAX  computer  facility  used.  The 
numerical  error  becomes  particularly  serious  as  the 
oscillating  frequency  becomes  large  where  the  extremely  thin 
boundary  layer  requires  more  grid  nodal  points  to  resolve 
the  flow  variables  in  the  vicinity  of  the  wall. 

Fig.  4-8  shows  that  as  the  Wormersley  number  gets 
large,  the  dimensionless  tidal  displacement  tends  to  the 
limit  of  2.0,  which  agrees  with  the  limit  of  1.0  in  the 
analytical  solution  given  by  Kurzweg  [18]  for  the  reason 
that  the  normalization  parameter  used  in  [18]  is  twice  as 
large  as  that  in  the  present  numerical  simulation. 

Fig.  4-9  shows  the  required  stimulating  axial 
pressure  gradient  used  in  the  present  study  for  a  pipe 
radius  R^  =  o.l  cm  and  water  (u  =  0.01  cm^/sec)  taken  as  the 
working  medium  versus  the  dimensional  tidal  displacement  ax 
in  cm  for  various  Wormersley  numbers  (namely,  oscillating 
frequency) .   It  is  evident  from  these  results  that  for  fixed 
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tidal  displacement,  the  required  axial  pressure  gradient  in 
the  large  Wormersley  number  case  is  orders  of  magnitude 
higher  than  that  in  the  small  Wormersley  number  case.  This 
may  eventually  put  some  constraint  on  the  use  of  very  high 
frequency  in  the  enhanced  heat  transfer  technique. 
Fortunately,  to  meet  the  tuning  condition,  the  required 
Wormersley  number  in  this  case  is  of  order  1.  As  already 
discussed  in  the  introduction,  in  order  to  gain  the  benefit 
of  axial  heat  transfer,  the  use  of  large  tidal  displacement 
is  always  preferred.  However,  such  an  increase  must  be 
limited  by  the  requirement  of  no  convective  net  mass 
transfer  occurring  between  two  reservoirs  and  may  also  be 
constrained  by  the  ability  of  the  device  to  withstand  the 
increase  in  the  exciting  pressure  gradient. 
Phase  Lags 

We  are  now  in  the  position  to  study  the  phase  lags  of 
the  Lagrangian  displacement  in  the  oscillatory  pipe  flow. 
Fig.  4-10  shows  the  phase  variations  (in  degrees)  along 
radius  for  Wormersley  number  varying  from  a  =  0.1  to  4. 
Some  numerical  results  are  also  shown  in  Tables  4-2  and  4-3. 
All  of  the  data  shown  in  Fig.  4-10  and  the  tables  have  the 
phase  angle  measured  relative  to  the  exciting  pressure 
gradient.  Two  features  can  be  seen  in  Fig.  4-10;  first,  in 
the  core  portion,  the  phase  lags  are  almost  equal  to  n/2 
when  the  Wormersley  number  is  small,  while  the  lags  are 
almost  n    when  the  Wormersley  number  is  large,  and  second. 
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the  phase  lags  vary  along  radius,  especially  in  the  boundary 
layer.  It  is  such  phase  lags  that  allow  the  existing 
temperature  gradient  in  the  very  thin  boundary  layer  of  the 
oscillating  pipe  flow  to  act  as  region  of  temporary  heat 
storage.   It  absorbs  heat  when  the  temperature  of  the  core 


Table  4-2 


Phase  Lags  Along  Radius 
(Working  Medium:  H2O,  Ax  =  10  cm) 


w^-- 


Nodal 

point 

K 


10 


a  =  0.1 


r?=r/Rl 


0 

.0000 

0 

.1667 

0 

.3333 

0 

.5000 

0 

.5833 

0. 

6667 

0. 

7500 

0. 

8333 

0. 

9167 

1. 

0000 

phase 


90.63 


90.63 


90.63 


90.63 


90.63 


90.63 


90.63 


90.63 


90.63 


0.00 


a    =    1. 


r?=r/Rl 


0 

.0000 

0 

.0366 

0 

.0957 

0 

.1913 

0 

.2593 

0 

3458 

0. 

4559 

0. 

5958 

0. 

7737 

1. 

0000 

phase 


102 

.06 

102 

.06 

102 

.06 

101 

.88 

101 

.70 

101. 

52 

101. 

16 

100. 

62 

99. 

72 

0.00 

a  =  2  • 


»?=r/Rl 


0 

.0000 

0 

.1095 

0 

.2164 

0 

.2836 

0 

.3620 

0 

4535 

0. 

5603 

0. 

6849 

0. 

8303 

1. 

0000 

phase 


134.28 


134.10 


133.38 


132.66 


131.76 


130.50 


128.70 


126.00 


122.58 


0.00 
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Table  4-2   —  continued 


Nodal 

point 

K 

a 

=  3 

a 

=  4 

a  = 

7 

»7=r/Rl 

phase 

r,=r/Rl 

phase 

r,=r/Rl 

phase 

1 

0.0000 

170.83 

0.0000 

179.82 

0.0000 

181.62 

2 

0.0917 

170.47 

0.1242 

179.28 

0.0970 

181.62 

3 

0.1818 

169.75 

0.2370 

178.02 

0.1892 

181.62 

4 

0.2703 

168.31 

0.3395 

176.04 

0.2769 

181.44 

5 

0.3573 

166.51 

0.4324 

173.70 

0.3604 

180.72 

6 

0.4427 

164.00 

0.5168 

171.00 

0.4397 

179.28 

7 

0.5267 

161.12 

0.5935 

168.12 

0.5152 

177.12 

8 

0.6091 

157.88 

0.6630 

165.24 

0.5869 

174.06 

9 

0.6901 

154.64 

0.7262 

162.18 

0.6552 

170.28 

10 

0.7696 

150.69 

0.7835 

159.12 

0.7201 

165.78 

11 

0.8478 

146.73 

0.8356 

156.06 

0.7818 

160.74 

12 

0.9246 

142.42 

0.8828 

153.18 

0.8406 

155.16 

13 

1.0000 

0.00 

0.9257 

150.30 

0.8964 

149.04 

14 

0.9646 

147.60 

0.9495 

142.56 

15 

1.0000 

0.00 

1.0000 

0.00 
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Table  4-2  —  continued 


Nodal 

Point 

K 


1^; 


a  =  10 


»7=r/Rl   phase 


1 

0.0000 

2 

0.0884 

3 

0.1740 

4 

0.2566 

5 

0.3364 

6 

0.4135 

7 

0.4880 

8 

0.5599 

9 

0.6295 

10 

0.6966 

11 

0.7616 

12 

0.8243 

13 

0.8849 

14 

0.9434 

15 

1.0000 

180 

.54 

180 

.72 

181 

.08 

181 

.80 

182 

.52 

182 

.70 

181. 

80 

179. 

64 

175. 

68 

170. 

46 

163. 

80 

156. 

06 

147. 

24 

137. 

52 

0.00 

a  =  14 


r)=r/Rl       phase 


0 

.0000 

0 

.0834 

0 

.1647 

0 

.2440 

0 

.3215 

0 

.3970 

0 

4709 

0, 

5428 

0. 

6130 

0. 

6815 

0. 

7484 

0. 

8136 

0. 

8773 

0. 

9394 

1. 

0000 

179.46 


179.10 


180.18 


184.50 


a  =  20 


r/=r/Rl   phase 


0.0000 


0.1635 


0.3032 


0.4227 


192.41 


189.53 


0.6120 


0.6865 


0.7502 


174.07   0.8046 


164.00 


153.20 


141.70 


130.19 


118.32 


0.00 


0.8512 


0.8909 


0.9249 


0.9540 


0.9788 


1.0000 


179 

.10 

180 

.54 

185 

.93 

182 

.34 

186 

.65 

184 

.14 

184. 

14 

181. 

98 

178. 

38 

173. 

71 

167. 

95 

161. 

12 

153. 

57 

146. 

01 

0,00 

■^ 
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Table  4-3   Comparison  of  Phase  Lags  With  Different 

Working  Mediums  (ax  =  10) 


Nodal 
Point 

K 

a    =    1. 

a  =  5. 

r,    =r/Rl 

Phase 

Lags 

r,    =r/Rl 

Phase 

Lags 

Water 

Mercury 

Water 

Mercury 

1 

0.0000 

102.06 

102.14 

0.0000 

181.62 

181.62 

2 

0.0161 

102.06 

102.14 

0.0714 

181.44 

181.26 

3 

0.0366 

102.06 

102.14 

0.1429 

180.90 

180.54 

4 

0.0626 

102.06 

102.14 

0.2143 

179.64 

179.46 

5 

0.0957 

102.06 

102.14 

0.2857 

178.02 

178.02 

6 

0.1378 

101.88 

102.14 

0.3571 

176.04 

175.86 

7 

0.1912 

101.88 

102.14 

0.4286 

173.34 

173.37 

8 

0.2593 

101.70 

102.14 

0.5000 

170.28 

170.11 

9 

0.3458 

101.52 

101.78 

0.5714 

166.50 

166.51 

10 

0.4559 

101.16 

101.42 

0.6429 

162.18 

162.20 

11 

0.5957 

100.62 

101.06 

0.7143 

157.50 

157.52 

12 

0.7737 

99.72 

99.98 

0.7857 

152.10 

152.13 

13 

0.0000 

0.00 

0.00 

0.8571 

146.16 

146.37 

14 

0.9286 

139.86 

139.90 

15 

1.0000 

0.00 

0.00 

W^ 
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slug  flow  is  higher  than  that  of  the  boundary  layer,  and  it 
releases  heat  as  the  core  temperature  is  relatively  lower. 
This  large  temperature  gradient  enhanced  by  the  existing 
velocity  phase  lags  allows  a  large  amount  heat  to  be 
conductively  transferred  radially  within  a  very  short  time 
and  subsequently  to  be  transferred  axially  by  a  convective 
coupling. 

Part  2.  The  Enhanced  Heat  Transfer  Investigation 
We  have  examined  some  mechanical  characteristics  in 
oscillating  pipe  flow,  and  compared  the  computed  solution  of 
the  velocity  field  with  Uchida's  solution.  The  results  for 
the  velocity  profiles  are  in  good  agreement.  Nevertheless, 
before  a  detailed  examination  of  the  thermal  field,  it  is 
first  necessary  to  test  the  current  developed  ETP  code  when 
applying  the  temperature  equations  (2-6) , (2-7) .  It  is  seen 
that  the  energy  equations  strongly  depend  on  the  velocity 
distribution  and  its  build-up  process,  so  that  one  can  use 
analytic  periodic  velocity  state  (Eqs.  1-4  and  1-5,  with  no 
build-up  process) ,  and  the  computed  velocity  (with  build-up 
process)  to  verify  the  correctness  of  the  resulting  thermal 
variables.  The  enhanced  heat  flux  is  a  function  of  both 
velocity  and  temperature  (Eqs.  2-34  and  2-35)  and  was  chosen 
for  a  comparison  of  the  analytic  and  numerical  results  of 
the  problem.  Part  A  of  table  4-4  shows  the  results  of  the 
computed  enhanced  axial  heat  flux  as  well  as  the  axial 
conduction  heat  flux  when  using  the  analytic  velocity 
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Table  4-4   The  Comparison  of  Enhanced  Heat  Flux  Using 
Numerical  Velocity  with  Heat  Flux  Using  Analytical  Velocity 
(Model  3,  Water-Glass,  Pr  =  7.03,  a  =  3) 


A.   Heat  Flux  Using  Analytic  Velocity 

AX 

(cm) 

(w/cm2°K) 

4>f 
(w/cm^-K) 

'f'V 
(w/cm2°K) 

(w/cm'*°K) 

0.9839 

0.0230 

0.0149 

0.0811 

0.0238 

2.9519 

0.1855 

0.0132 

0.0724 

0.0213 

4.9194 

0.3864 

0.0098 

0.0541 

0.0159 

6.8868 

0.5810 

0.0075 

0.0415 

0.0122 

7.8847 

0.6775 

0.0067 

0.0369 

0.0109 

9.8390 

0.8638 

0.0055 

0.0302 

0.0089 

B.   Heat  Flux  Using  Computed  Velocity 

1.0033 

0.0234 

0.0149 

0.0811 

0.0232 

3.0101 

0.1870 

0.0133 

0.0728 

0.0206 

5.0164 

0.3899 

0.0097 

0.0533 

0.0155 

7.0226 

0.5790 

0.0076 

0.0420 

0.0117 

8.0402 

0.6746 

0.0068 

0.0374 

0.0104 

10.0328 

0.8597 

0.0055 

0.0302 

0.0085 

*  Model  3,  Water-Glass 

*  (ft enhanced  axial  heat  flux 

*  <f>f axial  heat  flux  by  conduction  in  fluid 

*  <6v/ axial  heat  flux  by  conduction  in  pipe  wall 

*  yS  =  <^/AX^ 


87 
(Eq.  1-4)  for  a  =  1,  Pr  =  7.03,  pipe  radius  R^  =  0.1  cm  and 
glass  walls  of  thickness  AR  =  R2-  Ri  =  0.05  cm  (which  is 
almost  equal  to  the  boundary  layer  thickness  5  =  0.047  cm 
for  this  case)  with  Model  3.  There  were  22  points 
distributed  along  the  radius  (13  points  in  fluid),  and  the 
smallest  distance  next  to  the  wall  is  equal  to  0.166,  while 
along  the  axial  direction  101  point  were  used  (60  points  in 
the  central  pipe  section) .  Thirty  period  runs  were 
considered  in  order  to  insure  that  the  final  periodic  state 
was  achieved.  An  exception  was  the  Ax  =  1  cm  case  where 
only  23  period  runs  were  computed.  Under  the  exact  same 
conditions  described  above,  a  calculation  using  the 
numerical  velocity  with  build-up  process  was  tested.  Part  B 
of  table  4-4  gives  the  axial  heat  flux  results  using  the 
corresponding  numerical  generated  velocity  profile.  A 
comparison  of  results  shown  in  Table  4-4  is  quite  good.  One 
sees  that  the  disagreement  for  most  terms  is  less  than 
0.001.  Note  for  AX  of  10  cm,  the  ratio  (i>/4>f  is  about  15  and 
gets  larger  as  ax  is  increased  further. 
Periodic  Temperature  Build-up  in  Thermal  Pumping  Process 

The  purpose  of  this  facet  of  the  present  study  was  to 
examine  the  temperature  field  build-up  pattern  and  to 
determine  how  long  before  the  final  periodic  temperature 
state  can  be  reached.  The  author  also  intends  to  verify 
whether  or  not  there  indeed  exists  a  time  averaged  linear 
axial  temperature  distribution  within  the  thermal  pump  at 
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larger  time  as  assumed  in  existing  analytic  solutions  [15, 
16,  22]. 

Model  1  was  employed  to  study  the  temperature  build- 
up history  for  the  case  where  the  initial  velocity  was  zero 
everywhere  inside  the  pipe  as  well  as  at  the  boundary.  The 
initial  temperature  condition  was  assumed  to  be  identical 
with  the  cold  end  temperature  at  all  grid  points  except 
those  at  X  =  0  where  the  dimensionless  temperature  T  =  Tj^  = 
1.0 

Figs.  4-11,  4-12  and  4-13  show  the  computed 
temperature  build-up  history  at  Wormersley  number  a  =  1  at  x 
=  L/4,  and  L/2  (where  L  =  20  cm  is  the  pipe  length),  and 
either  on  the  axis  r  =  0  or  at  the  wall  r  =  R^.  The  tidal 
displacements  were  chosen  as  Ax  =  L/10,  L/4,  and  L/2, 
respectively.  In  these  figures,  the  ordinate  represents  the 
dimensionless  temperature,  while  the  abscissa  represents  the 
oscillatory  period  runs.  These  temperature  histories  show  a 
sinusoidal  variation.  Three  characteristics  can  be 
summarized:  1)  when  the  tidal  displacement  is  large,  the 
amplitude  of  the  temperature  variation  is  also  large,  and 
less  adjusting  time  is  needed;  2)  the  amplitude  of  the 
sinusoidal  variation  of  the  temperature  along  the  axis  is 
larger  than  that  near  the  wall  (see  Figs.  4-4,  4-5);  3) 
after  a  short  transient  process,  the  final  oscillatory  state 
is  reached.  The  time-averaged  temperatures  of  this 
oscillatory  state  at  a  given  x-position  have  the  same  value 


L 


f"*- 


•r'e 


t.eeBp  T 


89 


i.oee 


e.cee 


o.eee 


at  X  =  L/2,    r  =   0 


t.eeer- 


e.9ee 


e.eeo 


e 


at   X   =   L/2,    r  =  Ri 


J — I — L- 1 — I — \ — I — I I     I     ■     ■     I     I 


*7 


94 


141 


-J 1 1 1 L_L. 


188 


2J4 
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Fig  4-13    Temperature  Build-up  Process  In  Oscillating 

Flow  (Model  1,  q  =  1,  ax  =  10cm) 
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92 
for  all  tidal  displacements.  It  implies  that  the  final  time- 
averaged  temperature  depends  on  its  x-position  only.   For 
example,   the   final  periodic  dimensionless  time-averaged 
temperature  is  about  0.5  at  the  location  x  =  L/2,  while  it 
is  about  0.75  at  X  =  L/4  since  this  latter  position  is  much 
closer  to  the  hot  end.   From  this  observation,  one  may  well 
imagine  that  the  periodic  temperature   state  along  the 
longitudinal  direction  is  close  to  the  linear  function  form 
assumed  in  existing  analytic  studies  [17,  22].   This  feature 
will  be  further  demonstrated  later. 

Figs.  4-11,  4-12,  4-13  also  show  the  time  needed  for 
a  build-up  to  the  final  periodic  state  of  the  thermal  field. 
An  unpublished  conservative  analytical  estimation  of  the 
build-up  time  by  Kurzweg,  based  on  the  assumption  that  the 
oscillating  flow  satisfies  the  tuning  condition  is 


ti  =  l2/[wAx2  a  (a,Pr)  ] 


(4-1) 


Where   a  (a .  Pr)   =   3ceh/coAx2   is   a   function   of   thermal 
conductivity  and  Wormersley  number  [15,  16].   At  the  tuning 
point,   this  value  of  a   is  about  0.02.    if  the  tidal 
displacement  Ax  =  L/2  =  lo  cm  and  the  angular  frequency  .  = 
1  (for  water  with  R^  =  o.l  cm),  the  adjustment  time  t^  can 
be  estimated  from  equation  (4-1)  to  be  equal  to  200  sec,  or 
32  periods.    The  computed  results  show  that  the  actual 
adjusting  time  is  only  20  periods. 
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Fig  4-14    Temperature  Build-Up  Process  In  Steady  Flow 

(Model  2,  Uea  =  1.5  cm/sec^   ^ 


eq  =  1.5  cm/sec) 


94 


500 


AGO 


O 

•H 
Vi 

« 

0« 


300     • 


200 


/ 


/ 


/ 


/ 


7 


/ 


/ 
/ 

/ 

A Exciting 

'    Pressure  Gradient 


Build-up  Time 


!0 


15 


20 


50 


X 

^0  ^ 


30 


(U 


E 
20  <^ 

0) 

»^ 

.  3 
U 
U 
0) 


10 


25 


Tidal  Displacement  ax  (cm) 


Fig  4-15    Build-up  Time  Versus  Tidal  Displacement 

(Model  2,  Q  =  1) 


Ikj^ 


95 


Fig.  4-15  shows  the  build-up  time  versus  the  tidal 
displacement  and  applied  pressure  gradient  corresponding  to 
various  tidal  displacements  with  Model  2.   it  can  be  seen 
that  when  the  tidal  displacement  is  very  small,  the  build- 
up time  tends  towards  infinity,  so  that  the  so-called  final 
periodic  state  no  longer  exists.    The  final  periodic 
temperature  state  will  thus  be  established  only  by  heat 
conduction  and  has  a  complementary  error  function  erfc(x,t) 
=  1  -  erf(x,t)  temperature  build-up  history  [1,  H].   as  the 
tidal  displacement  gets  larger,  the  build-up  time  needed 
decreases,  and  when  ax  is  greater  than  about  L/4  (5  cm)  the 
build-up  time  is  almost  constant. 

In  order  to  compare  with  the  steady  Hagen-Poiseuille 
flow,  a  transient  temperature  history  of  such  a  flow  with 
Model  2  is  also  shown  in  Fig.  4-14.   Both  curves  represent 
the  temperature  computed  along  the  axis  (r  =  o) .   The  upper 
temperature  curve  is  at  x  =  L/4,  while  the  lower  one  is  at  x 
=  L/2.   It  is  clearly  seen  that  here  there  is  no  oscillatory 
pattern  at  all.   The  final  temperature  state  is  here  related 
to  the  Classical  Graetz  problem  dealt  with  in  great  detail 
in  the  heat  transfer  literature  [5,  H] . 
Temperature  nis-hr-jbutinn  in  ModPi  i 

The  final  periodic  temperature  distributions  for  the 
three  different  models  with  either  conducting  walls  or 
insulating  walls  have  been  obtained.  Fig.  4-16,  4-17  and  4- 
18  show  the  computed  periodic  final  temperature  distribution 
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Fig  4-16    Temperature  Distribution  in  Oscillating  Pipe 

Flow  (Model  1,  a  =  1,  AX  =  icm) 
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Fig  4-18   Temperature  Distribution  in  Oscillating  Pipe 

Flow  (Model  1,  a    =    1,    AX  =  5cin) 


99 
T(r,x,«)  in  oscillating  pipe  flow  with  Model  i.   The  marked 
numbers  1,  2,  3  represent  the  temperature  along  the  axis  (r 
=  0),  r  =  Ri/2,  and  at  r  =  R^  (wall),  respectively.   For 
these  computations,  the  Wormersley  number  was  taken  as  a  =  i 
and  the  tidal  displacement  was  chosen  as  Ax  =  i  cm,  2  cm  and 
5   cm,   respectively.     Neglecting   the   somewhat   erratic 
behavior  near  the  pipe  ends,  the  time-dependent  temperature 
shows  a  clear  linear  distribution  pattern  along  longitudinal 
direction  especially  for  the  relatively  small  tidal 
displacement  case,  of  say,  ax  =  1  cm.   if  one  focuses  on  one 
cross-section  (x  =    constant),  one  can  clearly  see  the 
temperature  near  the  boundary  and  that  at  the  axis  are 
alternately  higher  and  lower  relative  to  each  other  (Fig.  4- 
18)  within  a  period.   This  implies  that  a  very  large  radial 
conduction  heat  exchange  occurs  between  the  core  of  the  flow 
and  the  boundary  layer.   Note  that  the  radial  temperature 
gradient  becomes  very  large  during  parts  of  an  oscillating 
cycle. 

The  poor  behavior  of  the  temperature  distribution 
noted  near  both  ends  of  the  pipe  in  Fig.  4-18  at  tidal 
displacement  of  ax  =  5  cm  is  caused  by  numerical  dispersion 
which  is  a  common  shortcoming  of  the  second  order  numerical 
method  employed  [24,  28,  30,  37]  .  it  may  be  improved  by 
either  using  an  odd  order  numerical  method  such  as  the 
first-order  upwind  method,  flux  splitting  technique  [24]  or 
by  employing  alternative  models  which  are  able  to  avoid  any 
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discontinuity  of  the  physical  properties  at  both  ends,  such 
as  Model  3. 

Temperature  ni  stribiiti  on  in  Mnd^i  -> 

In  order  to  explore  some  alternative  possible  devices 
which  apply  the  enhanced  thermal  pumping  technique,  Model  2 

(Fig.  2-3)  was  developed  and  numerically  examined.  The 
temperature  distributions  in  Model  2  have  been  studied  and 
results  are  recorded  in  Figs.  4-19,  4-20,  4-21,  4-22  and  4- 
23  for  some  special  cases.  The  working  fluid  used  is  again 
water  and  similarly,  the  Wormersley  number  was  taken  as  a  = 

1.  The  tidal  displacements  used  were  Ax  =  i  cm,  5  cm,  10 
cm,  20  cm,  and  30  cm,  respectively.   The  marked  numbers  1, 

2,  3  in  those  figures  represent  the  temperatures  at  r  =  R^ 
(wall),  Ri/2  and  r=0  (i.e.,  along  the  axis).   it  can  be  seen 
that  the  temperatures  in  the  core  of  the  flow  and  that  in 
the  boundary  layer  in  these  figures  vary  periodically  in  the 
radial  direction  during  the  oscillations.   The  results  also 
show  a   linear  axial   temperature  variation  between  the 
heating  and  cooling  sources  when  the  tidal  displacement  is 
small  (AX  =  1  cm.  Fig.  4-19).   As  the  tidal  displacement 
becomes  larger  (ax  =  5  cm,  Fig.  4-20),  the  linear  axial 
temperature  variation  becomes  weaker.    The  difference  of 
temperature  in  the  core  of  the  flow  and  that  in  the  boundary 
layer  is  evident  and  shows  a  radially  alternating  pattern 
within  an  oscillation  period  just  as  seen  in  Model  1. 
However,  as  the  tidal  displacements  are  further  increased. 
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Fig  4-19    Temperature  Distribution  in  Oscillating  Pipe 

Flow  (Model  2,  q  =  1,  ax  =  Icm) 
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Fig  4-20   Temperature  Distribution  in  Oscillating  Pipe 

Flow  (Model  2,  a  =  1,  AX  =  5cm) 
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Fig  4-21        Temperature   Distribution   in  Oscillating  Pipe 

Flow    (Model   2,    a    =   1,    ax  =   10cm) 
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Fig  4-22    Temperature  Distribution  in  Oscillating  Pipe 

Flow  (Model  2,  a  =  1,  ax  =  20cm) 
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Fig   4-23        Temperature   Distribution    in   Oscillating   Pipe 

Flow    (Model    1,    q    =    1,    ax  =   30cin) 
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Fig  4-24    Temperature  Distribution  in  Steady  Flow 
(Model  2,  Uave  =  0.5-7.5  cm/sec) 


107 


as  shown  in  Fig.   4-21,  the  radial  temperature  gradient 
almost  no  longer  exists  except  directly  in  the  vicinity  of 
the  heat  source  and  sink  area.    it  is  obvious  that  the 
cross-section  averaged  temperature  in  the  pipe  tends  to 
become  constant  and  equal  to  the  mean  value  of  the  hot  and 
cold  source  region.   This  phenomenon  can  be  seen  even  more 
Clearly  as  the  tidal  displacement  is  further  increased 
(Figs.   4-22,   and  4-23).    One  can  conclude  from  such  a 
temperature  distribution  that  the  increase  of  axial  heat 
transfer  ability  with  increasing  tidal  displacement  is  most 
likely  to  be  limited.   At  the  same  time,  this  axial  heat 
transfer  capability  is  expected  to  increase  by  increasing 
either  the  individual  source  size  or  the  density  of  the 
alternatively  distributed  hot  and  cold  sources. 

This  unexpected  temperature  distribution  in  Model  2 
at  large  tidal  displacements  led  us  to  further  examine  the 
thermal   field   of   steady   pipe   flow   for   the   boundary 
conditions  used  in  Model   2.    The  computed  temperature 
distribution  with  an  equivalent  steady  averaged  velocity 
from  0.5  cm/sec  to  7.5  cm/sec,   corresponding  to  tidal 
displacement  from  1  cm  to  15  cm  at  an  oscillating  frequency 
Of  0,  =  i/sec  are  shown  in  Fig.   4-24.    The  equivalent 
velocity  is  established  from  the  relation  Ueq  =  ax<./2.   As 
could   have   been   anticipated,   large   radial   temperature 
gradients  in  this  flow  (similar  that  in  the  oscillating  flow 
with  large  tidal  displacement),   occur  only  at  the  heat 
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source  and  sink  areas,  while  the  radial  temperature  gradient 
in  the  rest  of  the  pipe  is  very  small.  The  constant 
temperature  difference  at  the  right-hand-side  differs  from 
that  on  the  left  hand  side  of  the  heat  source.  This 
temperature  difference  becomes  progressively  smaller  as  the 
equivalent  velocity  is  increased  (i.e.,  the  fluid  heating 
becomes  less) . 
Temperature  Distribution  in  Model  3 

Model  3  was  developed  to  examine  the  performance  of 
the  thermal  pump  without  convective  heat  exchange  at  the 
ends  (in  the  chambers)  .    That  is,  the  heat  is  added  and 
withdrawn  by  conduction  through  the  pipe  end  section  (see 
Fig.  2-4).   The  computed  temperature  distribution  in  this 
model  is  shown  in  Fig.  4-25,  and  a  magnified  view  of  the 
temperature  in  the  central  section  is  shown  in  Fig.  4-26. 
The  Wormersley  number  used  is  again  a  =  1  and  the  tidal 
displacement  Ax  =  10  cm.   Two  features  are  observed.   First, 
in  contrast  to  the  temperature  computed  in  Model  1,  the 
temperature  obtained  here  shows  a  very  smooth  pattern  at  the 
central  section  ends  (i.e.,  at  x  =  5L  and  6L  in  Fig.  2-4). 
Second,  the  linear  axial  temperature  distribution  with  time- 
dependent  radial  temperature  gradient  alternating  in  sign  in 
the  central  portion  of  the  pipe  (namely,  5L  <  x  <  6L,  in 
Fig.  2-4)  is  retained  even  when  the  tidal  displacement  is  as 
large  as  ax  =  10  cm  (Fig.  4-26) . 
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Fig  4-25   Temperature  Distribution  in  Oscillating  Pipe 
Flow  (Model  3,  a  =  1,  AX  =  lOcm,  L=20cin) 
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Fig  4-26   Magnified  View  Of  Temperature 
in  the  Central  Pipe  Section 
(Model  3,  Q  =  1,  AX  =  10cm) 


Ill 

Heat  Flux  Versus  Tidal  Displacement  fModel  2) 

As  mentioned  in  the  previous  section,  a  limit  on  the 
axial  heat  transfer  in  oscillating  flow  with  Model  2  exists, 
and  the  heat  axially  transferred  in  the  steady  Hagen- 
Poiseuille  flow  with  the  same  model  could  be  expected  as  the 
upper  limit  of  the  corresponding  oscillating  flow. 

This  has  been  verified  by  a  numerical  calculation 
with  water  as  the  working  medium  and  the  results  are 
recorded  in  Table  4-5,  and  4-6,  and  also  plotted  in  Fig.  4- 
27.   The  Wormersley  number  used  in  the  current  test  is  a  = 
1,  and  the  tidal  displacement  for  the  oscillating  flow 
varied  from  l  cm  to  30  cm.   The  corresponding  equivalent 
average  velocity  varied  from  0.5  cm/sec  to  15  cm/sec  for  the 
steady  flow  case.   In  Fig.  4-27  the  ordinate  represents  the 
heat  flux  released  from  the  hot  source,  the  abscissa  is  the 
tidal  displacement  or  the  corresponding  equivalent  velocity 
for  steady  flow.    The  solid  curve  shows  the  heat  flux 
computed  for  the  oscillating  flow,  while  the  dashed  curve 
shows  the  axial  heat  flux  in  steady  flow.   It  can  be  seen 
that  the  slope  of  the  heat  flux  curve  in  the  oscillating 
pipe  flow  case  is  smaller  than  that  in  steady  pipe  flow  in 
general  and  particularly  when  the  tidal  displacement  is 
small  (i.e.,  ax  <  3  cm).   The  flux  increases  rapidly  as  the 
tidal  displacement  gets  larger  (ax  =  4-6  cm) .   Still  further 
increasing  the  tidal  displacement  does  not  lead  to  a  further 
large  increase  of  the  axial  heat  flux.   In  fact  the  slope 
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begins  to  decrease  after  the  tidal  displacement  becomes 
larger  than  about  Ax  =  8  cm.  Eventually  this  slope  tends  to 
zero  as  ax  increases  still  further.   For  the  steady  flow. 


Table  4-5  Enhanced  Heat  Flux  via  Tidal  Displacement 
(Model  2,  oscillating  flow  Pr  =  7.03,  a   =  3) 


AX 

(cm) 

Periods  to 

final 
oscillation 

state  * 

1.000 

75 

1.962 

15 

2.943 

7 

3.924 

6 

4.905 

5 

6.867 

5 

7.848 

5 

9.843 

5 

15.233 

5 

19.920 

5 

30.100 

5 

Max.  Pressure  gradient 
h=l/p\dp/dx\     (cm/s2) 


4.86 


8.372 


12.588 


16.744 


20.930 


29.300 


33.488 


42.000 


65.000 


85.000 


126.000 


heat  flux  <f) 
(w/cm2 - K) 


0 

.03 

0 

.136 

0 

.345 

0 

,683 

1 

.051 

1 

.377 

1. 

420 

1. 

484 

1. 

545 

1. 

522 

1. 

530 

*  build-up  time  can  be  estimated  by  equation  (4-1) 
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Table  4-6  Enhanced  Heat  Flux  in  Steady  Flow 


Ax 

(cm) 

Max. 
Pressure  gradient 
A=i/p  ap/ax| 
(cm/s^) 

Ueq=0.5wAx 
(cm/s) 

heat  flux  4> 
(w/cm^-K) 

0.500 

2.43 

0.25 

0.462 

1.000 

4.86 

0.50 

1.354 

2.500 

10.50 

1.25 

1.567 

3.000 

12.56 

1.50 

1.631 

5.000 

21.00 

2.50 

1.775 

8.000 

33.50 

4.00 

1.664 

10.000 

42.00 

5.00 

1.786 

15.000 

65.00 

7.50 

1.779 

the  slope  of  the  heat  flux  curve  is  quite  large  and  almost 
constant  at  all  equivalent  velocities  less  than  0.5  cm/sec, 
but  it  rapidly  decreases  when  the  equivalent  velocity 
becomes  larger.  Similar  to  the  oscillating  flow  case,  the 
heat  flux  ^  as  the  velocity  increases  reaches  a  limit.  This 
limit  is  about  20  percent  greater  than  that  for  the 
oscillatory  flow. 

The  reason  of  the  slope  variation  in  the  oscillating 
flow   can   be   explained   as   follows:   At   small   tidal 
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displacement,  k^       is   quite   small,   and   the   molecular 
conduction  is  dominant,  so  that  only  a  very  small  quantity 
of  heat  will  be  transferred  axially.   However,  as  the  tidal 
displacement  increases.   Kg  =  Au,Ax2  becomes  considerably 
larger   and   the   enhanced   heat   transfer   mechanism   now 
dominates  axial  conduction.   The  limit  on  axial  heat  flux 
observed  at  large  ax  is  due  to  a  bottleneck  which  develops 
at  the  heat  supply  and  removal  sections  due  to  insufficient 
conduction  heat  transfer. 

In  the  steady  flow,  the  fluid  particles  are  always 
fresh  ones  with  no  residual  heat  in  them  so  the  radial 
temperature  gradient  is  relatively  larger  than  that  in  the 
oscillating  pipe  flow  case.  The  larger  radial  temperature 
gradient  increases  the  heat  supply  and  removal  capability  so 
as  to  ensure  an  axial  flux  limit  in  excess  of  that  possible 
for  oscillating  pipe  flow  in  this  model. 
Influence  of  The  TbprmodvnaTnin  Pr-pperti^^  (m.,h^i  .^ 

It  is  well  known  that  the  thermodynamic  properties  of 
the  working  fluid,  especially  the  viscosity,  are  functions 
of  the  local  temperature.  Some  thermodynamic  properties  of 
water  for  the  temperature  range  of  0=C  to  200»C  are  listed 
in  Table  4-7.  One  can  see  that  even  over  such  a  small 
temperature  range  the  thermodynamic  properties  of  water  vary 
considerably.  For  instance,  at  T  =  o  =  C,  the  kinematic 
viscosity  .  =  1.788  lO'S  M2/Sec,  and  the  corresponding 
Prandtl  Number  is  Pr  =  13.6,  while  at  temperature  T  =  ioo° , 
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Table  4-7   The  Influence  of  Properties  of 
Water  on  the  Enhanced  Axial  Heat  Flux 


T 
C 


Pr 


200 

0.937 

160 

1.099 

100 

1.740 

80 

2.22 

60 

3.02 

40 

4.34 

20 

7.02 

0 

13.60 

g/citi- 


0 

.8668 

0 

.9097 

0 

.9606 

0 

.9741 

0, 

9855 

0, 

9946 

1. 

0000 

1. 

0020 

ws/g-K 


4 

.505 

4 

.417 

4 

.216 

4 

.196 

4. 

184 

4. 

178 

4. 

182 

4. 

218 

u    10~2 
cm^/sec 


0 

.160 

0 

.173 

0 

.294 

0 

.364 

0. 

478 

0. 

658 

1. 

006 

1. 

788 

K 

w/cm" K 


0.665 


0.680 


0.680 


0.668 


AX 


cm 


0.628 


0.597 


0.552 


10 

.031 

11 

.499 

10 

.029 

10 

.030 

10. 

030 

10. 

041 

9. 

374 

10. 

030 

4> 
at  a=l 
w/cm^'K 


1 

.025 

1 

.249 

1 

.500 

1 

.525 

1. 

706 

1. 

780 

1. 

476 

0. 

747 

}i 


*   Properties  are  in  the  saturated  state  [11]. 

n  =  0.294  10-6  M2/sec  and  Pr  =  1.74.   It  varies  almost  an 
order  of  magnitude  over  this  range. 

The  present  study  did  not  intend  to  investigate  the 
influence  of  property  changes  in  full  detail;  however,  the 
performance  of  axial  heat  flux  versus  water  properties  at 
temperature  from  0°C  to  200° C  has  been  examined. 

Though  Model  2  may  not  be  an  ideal  device  (at  large 
AX)  for  the  application  of  enhanced  thermal  pumping,  it  can, 
however,  be  a  good  model  for  demonstrating  the  influence  of 
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the  variation  of  the  fluid  properties  on  the  heat  transfer 
near  the  source  area  in  oscillating  pipe  flow.  We  used  the 
tidal  displacement  of  ax  =  10  cm,  which  according  to  Fig.  4- 
27,  lies  in  the  maximum  axial  heat  flux  range.  The  enhanced 
heat  flux  of  the  final  periodic  state  is  shown  in  Fig.  4-28 
by  the  solid  curve  as  a  function  of  the  Prandtl  number  for 
water  at  different  temperatures  and  also  in  the  last  column 
of  Table  4-7.  The  dashed  curve  shown  in  Fig.  4-28  is  the 
time  needed  to  build-up  to  the  final  oscillation  state.  The 
Wormersley  number  used  in  the  calculations  was  a  =  1. 

According  to  the  definition  of  Wormersley  number 
given  by  equation  (1-10),  for  a  fixed  a  and  pipe  radius  R^, 
the  angular  velocity  w  is  directly  proportional  to  the 
kinematic  viscosity  «/ .    This  also  implies  that  the  fluid 
boundary  layer  thickness  is  constant  (Eq.  1-2) .   So,  the 
solid  curve  in  Fig.   4-28,   in  fact,  mainly  reflects  a 
relationship   between   the   enhanced   heat   flux   and   the 
oscillating  frequency  «.    It  is  found  that  a  maximum 
enhanced   heat   flux   can  be   obtained   if   the  water   of 
approximate   temperature   T   =   40°  C   is   employed.     The 
corresponding  angular  velocity  at  this  temperature  level  is 
w  =  0.658  radian/sec.   Physically,  the  existence  of  a  peak 
value  in  Fig.  4-28,  may  be  explained  as  follows:  At  higher 
temperature,  say,  T  =  100°  C,  the  kinematic  viscosity  v    is 
quite  small  and  hence  the  angular  velocity  w  is  too  slow  to 
transport   large   amounts   of     heat   axially.     As   the 
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temperature  decreases,  i/  gets  larger,  and  so  does  the 
angular  velocity  w.  Thus  more  heat  can  be  expected  to  be 
transported  and  this  eventually  reaches  the  limit  mentioned 
in  the  last  section.  It  seems  that  a  similar  flattening 
pattern  in  the  heat  flux  should  also  appear  if  the 
temperature  is  further  decreased  and  hence  the  frequency  w 
further  increased.  However,  as  temperature  decreases 
further  the  thermal  conductivity  declines  (Table  4-7)  so 
that  it  weakens  the  heat  supply  and  removal  capability  at 
the  sources  and  hence  makes  the  axial  heat  flux  curve 
decline  from  its  maximum. 
Heat  Flux  Versus  Tidal  Displacement  (Model  1) 

As  indicated  by  Kurzweg  in  earlier  analytic  studies 
[16,  19,  20],  the  enhanced  axial  heat  transfer  versus  the 
oscillating  flow  is  directly  proportional  to  the  square  of 
the  tidal  displacement.  A  series  of  numerical  tests  with 
Model  1  were  run  to  check  this  prediction  and  to  study  the 
influence  of  various  fluid  properties  and  boundary 
conditions  on  this  Ax^  effect. 

The  computations  were  separated  into  two  groups:  in 
the  first  group,  the  Wormersley  number  was  chosen  very  close 
to  the  tuning  point,  i.e.,  a  =  1.  Water  was  chosen  as  the 
working  fluid,  and  both  insulating  and  conducting  wall 
conditions  (glass  and  steel)  were  considered.  The 
computational  results  are  shown  in  Fig.  4-29.  In  the  second 
group  of  numerical  studies,  the  Wormersley  number  was  chosen 
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Fig  .-29   Heat  Flux  versus  Tidal  Displace„ent 

(Model  1,    a    -    I,    pr   -    7.03) 
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Fig  4-30    Heat  Flux  versus  Tidal  Displacement 

(Model  1,  a    =  3,  Pr  =  7.03) 
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as  a  =  3,  which  deviates  from  the  tuning  point.  Again  water 
was  chosen  as  the  working  medium.  Similarly,  both 
insulating  wall  and  conducting  (glass  only)  wall  cases  were 
investigated.  The  results  are  shown  in  Fig.  4-30.  The 
solid  curve  (Fig.  4-29)  shows  the  enhanced  axial  heat  flux 
in  oscillating  pipe  flow  versus  tidal  displacement  for  the 
insulating  wall  case  and  the  dashed  curve  represents  the 
enhanced  axial  heat  flux  obtained  with  a  steel  wall.   The  " 

.  .  .  n    curve  shows  the  enhanced  axial  heat  flux 

computed  for  a  glass  wall.    These  three  curves  clearly 
confirm  the  quadratic  behavior, 

^  =  k^«Ax2 

where  k^  is  a  constant.   For  this  specific  case,  it  is  found 
that 

k^  =  0.0056         for  the  insulating  wall  case, 
k^  =  0.0086         for  the  steel  wall  case, 
k^  =  0.0098         for  the  glass  wall  case. 
The  computational  results  for  the  Wormersley  number  a 
=  3  case  are  shown  in  Fig.  4-30.   There  the  solid  curve 
represents  the  enhanced  axial  heat  flux  in  oscillating  pipe 
flow  with  insulating  wall,  the  dashed  curve  shows  the  glass 
wall   results.     Similarly,   the   coefficient   k^   in   the 
quadratic  formula  is  found  to  have  the  value 

k^  =  0.0233         for  the  insulating  wall  case, 
k^  =  0.0317         for  the  water-glass  combination. 
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It  is  obvious  that  the  existence  of  conducting  walls 
tends  to  increase  heat  storage  capacity  and  hence  enhance 
the  radial  and  axial  heat  transfer  process.   Thus  one  can 
expect  a  larger  heat  flux  for  conducting  walls  than  for 
insulating   walls   at   the   same   tidal   displacement   and 
Wormersley  number.  It  can  be  seen  from  Fig.  4-29  and  Fig.  4- 
30  that  higher  Wormersley  number  (a  =  3)  will  generally  lead 
to  a  larger  axial  heat  flux.   In  these  specific  cases  higher 
wormersley   number   simply   implies   higher   oscillating 
frequency  a,.  However,  it  does  not  mean  that  the  tuning  point 
shifts  to  a  Wormersley  number  equal  to  3.   Because  here  the 
enhanced  axial  heat  flux  is  a  dimensional  quantity,  if  this 
quantity  is  divided  by  the  oscillating  frequency,  one  can 
clearly  see  that  this  new  quotient  is  larger  in  the  a  =  i 


case. 


If  one  compares  the  computational  results  obtained  by 
using  Model  1  in  this  section  with  the  previously  obtained 
results  for  Model  2  as  shown  in  Fig.  4-27,  one  may  at  first 
glance  feel  somewhat  puzzled  at  discrepancies  between  the 
results.     Note,   however,   that   there   is   an   essential 
distinction  between  Model  i  and  Model  2.   in  Model  1,  the 
fluid  elements  have  infinite  heat  exchange  capability  at 
both  ends,  in  fact,  all  the  new  fluid  elements  with  constant 
temperature  (i.e.,  the  reservoir's  temperature)  enter  the 
pipe  at  both  ends  during  every  oscillating  period.   This 
implies  that  direct  convective  heat  exchange  occurs  at  these 
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ends  during  each  oscillation  so  as  to  match  the  fast 
enhanced  axial  heat  transfer  possible  within  the  pipe.   On 
the  other  hand,  for  the  Model  2  case,  there  is  no  such 
infinite  heat  exchange  capability  at  the  heat  source  and  the 
sink  sections.   in  fact,  the  heat  being  supplied  and  removed 
is  purely  by  molecular  heat  conduction  at  the  interface 
between  of  the  thermal  sources  and  the  fluid  boundary  layer. 
This   conductive   heat   supply   and   removal   mode   is   not 
compatible  with  a  very  fast  enhanced  axial  heat  transfer 
within  the  connecting  pipe  flow  and  actually  prevents  the 
thermal  process  from  working  at  peak  efficiency.    This 
comparison  suggests  to  us  that  one  should  consider  a  strong 
and  preferably  turbulent  convective  motion  in  both  end 
reservoirs  as  a  means  to  increase  the  efficiency  of  the 
enhanced  thermal  pumping  technique. 
The  Influence  of  Wall  Thickn^Qg 

As  mentioned  above,  the  conducting  wall  will,  in 
general,  make  the  enhanced  axial  heat  transfer  process  more 
efficient  since  the  existence  of  a  conducting  wall  creates 
an  additional  heat  storage  capability  near  the  fluid-solid 
interface.  This  additional  heat  storage  supplements  the 
storage  capacity  of  the  boundary  layer  and  hence  strengthens 
the  enhanced  heat  exchange  process  of  the  whole  system. 

The  largest  radial  distance  which  the  heat  can  be 
transferred  into  the  wall  is  the  wall  thermal  penetrating 
thickness  s^.       This  thickness  is  a  function  of  oscillating 
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frequency  as  well  as  the  thermal  properties  of  the  solid 
wall.  As  mentioned  in  the  introduction  chapter  (Fig.  2-1), 
the  future  engineering  applications  of  the  thermal  pumping 
technique  probably  will  require  a  bundle  of  capillary  tubes 
connecting  the  cold  and  hot  fluid  reservoirs.  in  view  of 
space  limitations  and  economics,  one  wants  to  minimize  the 
capillary  wall  thickness  used  in  the  thermal  pump  yet  still 
be  able  to  make  full  use  of  the  wall  heat  storage  capacity. 

Using  Model  1,  a  numerical  investigation  of  enhanced 
axial  heat  flux  in  oscillating  pipe  flow  in  a  water-glass 
combination  for  various  wall  thicknesses  has  been  performed 
and  the  computed  results  recorded  in  Fig.  4-31.    m  this 
study,  the  tidal  displacement  was  taken  as  ax  =  5  cm  (  = 
L/4),  the  Wormersley  number  a  =  i,  and  inner  pipe  radius  R^ 
=  0.1  cm.   water  at  20°  c  was  used  as  the  working  medium. 
If  we  define  the  effective  heat  flux  as  the  enhanced  axial 
heat  flux  computed  not  only  over  the  cross-section  of  the 
fluid  in  the  pipe  but  also  including  the  solid  wall  area,  we 
can  see  from  Fig.  4-31  that  as  wall  thickness  increases,  the 
effective  heat  flux  curve  first  increases  with  increasing 
wall   thickness   R2   -   r,   b^t   soon   starts   to   decline. 
Apparently  there  is  an  optimum  value  of  the  wall  cladding 
thickness.   For  this  specific  case,  it  has  been  found  that 
this  optimum  wall  thickness  is  about  10  percent  of  the  pipe 
diameter,   it  also  implies  that  the  penetrating  thickness  s^ 
into  the  wall  in  this  case  is  about  0.02  cm,  which  is  20 
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Fig  4-31   Influence  of  Wall  Thickness  on  Axial  Heat  flux 

(Model  1,  Water-Glass,  a  =  1,  ax  =  Scm) 
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percent  of  the  pipe  inner  radius.    The  decline  of  the 
effective  heat  flux  curve  is  due  to  the  "unpenet rated"  outer 
part  of  the  wall  which  plays  a  useless  role  in  the  enhanced 
heat  axial  transfer  process. 
The  Influence  of  Pipe  ni  stti^i-ot- 

Another  very  important  aspect  of  the  heat  transfer 
technique  under  consideration  is  the  performance  of  heat 
flux  versus  pipe  diameter.  A  numerical  simulation  with' 
Model  3  has  been  carried  out  to  examine  the  pipe  diameter 
effects.  For  the  test  the  Wormersley  number  was  taken  as  a 
=  3,  the  tidal  displacement  ax  =  lo  cm,  and  the  fixed  wall 
thickness  equal  to  0.05  cm.  Once  again  20»C  water  was  taken 
as  the  working  fluid  and  the  wall  material  was  glass. 

From  the  definition  of  Wormersley  number   (1-10), 
namely,  a     =     ^^Jl^ ,      it  is  clear  that  the  oscillating 
frequency  is  inversely  proportional  to  the  second  power  of 
pipe  inner  radius  for  fixed  Wormersley  number  and  directly 
proportional  to  the  kinematic  viscosity  of  the  working 
fluid.   on  the  other  hand,  the  enhanced  axial  heat  flux  is 
an   increasing   monotonic   function   of   the   oscillating 
frequency.    This  implies  that  an  increase  in  inner  pipe 
diameter  will   surely  lead  to  a  decline  of  oscillating 
frequency  and  hence  the  enhanced  axial  heat  flux,  if  both 
Wormersley  number  and  the  kinematic  viscosity  are  held 
constant.   The  numerical  results  are  shown  in  Fig.  4-32. 


-I 


128 


2. Or 


Pipe  diameter  R^   (cm) 


Fig  4-32   Influence  of  Pipe  Diameter  on  Heat  flux  for  Fixed 
Frequency  (Model  3,  Water-glass,  a  =  3,  Ax  =  10cm) 
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Fig  4-33  Typical  Iso-Temperature  Contour 

in  Oscillating  Pipe  Flow 
(Model  3,  a  =  3,  AX  =  locm,  Ri=0.1cin,  R2=0.5cm) 
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It  will  be  noted  that  the  4>  versus  R^  curve  does  not 
show  a  decrease  in  enhanced  axial  heat  flux  as  the  radius 
becomes  very  small.  It  is  expected  from  the  early  analytic 
work  that  there  exists  an  upper  limit  of  the  enhanced  axial 
heat  flux  as  the  inner  pipe  diameter  becomes  extremely 
small.  This  is  because  that  when  the  pipe  cross-section 
area  is  further  decreased  such  that  the  central  slug  flow 
volume  which  transports  heat  back  and  forth  axial ly  becomes 
very  small,  the  heat  transfer  mechanism  of  the  thermal 
pumping  technique  is  destroyed  since  large  radial 
temperature  gradients  are  no  longer  possible. 

Fig.  4-33  shows  iso-temperature  contours  for  a  pipe 
diameter  equal  to  0.5  cm.  The  iso-temperature  contours  show 
a  periodic  temperature  pattern  within  an  oscillation  period. 
It  is  noted  that  the  contours  at  phase  angles  wt  =  0°  and 
180°  or  90°  and  270°  are  not  precisely  anti-symmetric. 
Variation  of  Axial  Temperature  Gradient  in  Model  3 

In  order  to  overcome  the  numerical  instabilities 
observed  near  the  pipe  ends  in  Model  1  and  to  examine  the 
effect  of  no  convective  heat  exchange  at  the  pipe  ends, 
Model  3  which  is  shown  in  Fig.  2-4,  was  developed  and 
studied. 

The  first  numerical  study  carried  out  with  this  later 
model  was  to  examine  the  variation  of  the  time-averaged 
axial  temperature  gradient  within  the  central  pipe  section 
as  the  tidal  displacement  was  increased.   This  time-averaged 
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Fig  4-34   Variation  of  Temperature  Ti  and  T2  versus  AX 

(Model  3,  Water-Glass,  a  =  3) 
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axial  temperature  gradient  variation  can  be  calculated  by 
measuring  the  time-averaged  temperatures  at  the  joint  points 
(X  =  5L,  and  6L,  see  Fig.  2-4).   The  computed  results  are 
shown  in  Fig.  4-34.    m  the  calculation,  the  Wormersley 
number  was  chosen  as  «  =  i  and  the  working  medium  is  water 
within  a  glass  pipe.   m  the  figure,  T2  is  the  time-averaged 
temperature  computed  at  the  left  joint  (x  =  5L)  which  is 
close  to  the  hot  source,  while  Ti,  at  right  joint,  is  close  ■ 
to  the  cold  source.   it  is  clearly  seen  that  T2  decreases, 
while  Ti  increases  as  the  tidal  displacements  increase.   It 
seems  that  both  T^  and  T2  tend  to  the  dimensionless  mean 
value   (i.e.,   0.5)   if  the  tidal  displacement  is  greatly 
increased.    This  implies  that  the  time-averaged  axial 
temperature  gradient  in  the  central  section  will  decrease 
and  that  the  enhanced  axial  heat  flux  will  be  weakened 


Table  4-8   Variation  of  the  Axial  Temperature  Gradient 
versus  Wormersley  Number  (Water-GlaL,  Ax  2  ?o  cm? 
Ui,  T2 are  dimensionless  TemDeratuT-*.c:^ 


Temperatures) 
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although  still  increasing  with  increasing  Ax.   However,  the 

numerical  results  show  that  such  a  weakening  of  the  axial 

temperature  gradient  can  be  improved  at  higher  oscillating 

frequency  (Table  4-8).   This  is  because  at  higher  a,  the  pipe 

flow  develops  a  much  thinner  boundary  layer  in  the  vicinity 

of  the  wall  and  hence  a  larger  radial  temperature  gradient 

which  increases  the  capability  of  heat  supply  and  removal 

from  the  extended  sections.   In  fact,  the  numerical  results 

also  show  that  only  very  narrow  areas  of  the  extended 

section  are  involved  in  supplying  and  removing  heat  when  the 

oscillating  frequency  is  large  (Fig.  4-25) . 

Comparison  of  Enhanced  Osm' 1  i.atorv  Hp;^.^  Tr-^r.of^-r 
and  Heat  Conduction  "• 

Fig.  4-35  shows  the  enhanced  axial  heat  flux  computed 
in  the  central  pipe  section  in  Model  3  and  the  axial 
molecular  conduction  heat  flux  both  in  the  fluid  and  in  the 
wall  versus  tidal  displacement.   The  Wormersley  number  in 
the  test  case  is  chosen  as  a  =  3,  the  working  medium  is 
.water,  and  the  solid  wall  material  is  glass.    curve  ^ 
represents  the  variation  of  enhanced  axial  heat  flux  versus 
the  tidal  displacement.   This  is  seen  to  be  smaller  than  the 
heat  flux  due  to  pure  molecular  conduction  either  in  the 
wall  (curve  <f>^)     or  in  the  fluid  (curve  ^f)  when  the  tidal 
displacement  is  very  small  (ax  less  than  1  cm  in  this  case) /  ' 
however,    ^    increases   dramatically   when   the   tidal 
displacement  becomes  larger.   m  contrast,  the  pure  axial 
molecular  conduction  declines   slightly  due  a  weakening  of 
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Enhanced  Axial  Heat  Flux 
Axial  Heat  Conduction  Rate  in  Wall 
Axial  Heat  Conduction  Rate  in  Fluid 
^/Ax-^  (Watt/°K) 


idal  Displacement  ax  (cm) 


Fig  4-35 


Comparison  of  Enhanced  Heat  Transfer  and  Ro.^ 
Conduction  in  Oscillating  P?pe  Flow 
(Model  3,  Water-Glass,  a  iPJ^^^^ 
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the  axial  temperature  gradient  in  the  central  pipe  section 
at  large  Ax  as  discussed  in  the  last  section. 

The  ratio  of  the  enhanced  heat  flux  to  the  square  of 
the  tidal  displacement  is  also  shown  in  the  same  figure 
(curve  fi).  It  indicates  a  declining  value  of  p  as  the  tidal 
displacement  increases.  This  indicates  that  the  2nd  power 
relationship  between  the  enhanced  axial  heat  flux  4>  and 
tidal  displacement  ax  has  also  been  weakened  (i.e.,  lower 
k^)  if  large  tidal  displacements  are  used. 
Enhanced  Heat  F1uy  as  a  Fnnr.tinn  pf  worT^^r-^iev  Number- 

To  study  the  influence  of  oscillating  frequency  on 
the  enhanced  axial  heat  transfer,  a  series  of  numerical 
tests  were  performed  with  Model  3  for  various  combinations 
of  water-steel  and  mercury-steel  cases.   Fig.  4-36  shows  the 
enhanced  axial  heat  flux  and  the  heat  flux  by  pure  molecular 
heat  conduction  both  in  the  fluid  and  in  the  wall  as  a 
function   of   Wormersley   number   (i.e.,   a   function   of 
oscillating  frequency  .,  for  fixed  pipe  radius  and  working 
fluid) . 

The  solid  curves  show  the  computed  results  in  the 
water-glass  wall  case,  while  the  dashed  curves  represent  the 
numerical  solution  obtained  with  the  mercury-steel  wall 
combination.  The  tidal  displacement  ax  used  in  these 
calculations  is  fixed  and  equals  10  cm.  The  subscripts  have 
the  same  meaning  as  those  used  in  Fig.  4-35. 
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Fig  4-36   Variation  of  Axial  Heat  Flux  Versus  Wormersley 

Number  (Model  3,  Ax  =  10cm) 
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It  can  be  seen  from  Fig.  4-36  that  both  the  enhanced 
axial  heat  flux  4>    and  axial  heat  flux  due  to  pure  molecular 
conduction  in  either  the  wall  (^^)  or  in  the  fluid  (^f) 
increase  as  the  oscillating  frequency  gets  larger.   It  is 
assumed  that  this  is  partially  due  to  the  recovery  of  the  \ 
axial  temperature  gradient  when  the  oscillating  frequency 
increases  as  discussed  in  the  previous  section.   it  can  been 
also  seen  that  the  enhanced  axial  heat  flux  is  larger  than 
the  axial  heat  flow  by  pure  molecular  conduction  in  either 
case.   The  difference  between  the  enhanced  axial  heat  flux 
and  the  axial  heat  flux  due  to  pure  molecular  conduction  in 
the  water-steel  wall  case  is  smaller  than  that  in  the 
mercury-steel  case.    This  implies  that  the  water-steel 
combination  does  not  take  optimum  advantage  of  the  enhanced 
thermal  pumping  process.   The  axial  heat  flux  due  to  pure 
molecular  heat  conduction  in  the  fluid,  in  general  is  very 
much  smaller  than  the  enhanced  axial  heat  flux  and  hence  is 
indeed  negligible  as  assumed  in  the  theoretical  analysis  of 
this  problem. 

Note  that  if  the  dimensional  enhanced  heat  flux  is 
divided  by  the  corresponding  oscillating  frequency,  we  will 
obtain  the  tuning  curves  which  will  be  discussed  below. 
Tuning  Curves 

In  Fig.  4-37,  the  left  curve  shows  the  variation  of 
the  ratio  of  the  enhanced  axial  heat  flux  <f>  to  the 
oscillating  frequency  .  in  the  water-glass  wall  case,  while 
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Fig  4-37    Computed  Tuning  Curves 
(Model  3,  H20-Steel  and  Hg-steel,  Ax  =  lOcm) 
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the  right  curve  represents  the  mercury-glass  wall  case  as 
marked  in  the  diagram.  It  can  be  seen  that  both  curves  show 
the  existence  of  the  so-called  tuning  effect.  The  optimum 
value  of  R^  for  the  water-glass  case  (Pr  =  7.03,  w  =  1/sec) 
is  found  equal  to  0.1  cm,  and  the  corresponding  Wormersley 
number  is  about  1. 

Fig.  4-38  shows  some  analytical  results  obtained  in 
reference  [20,  21],  There  the  optimum  value  of  a  for  a 
Prandtl  number  Pr  =  10  is  close  to  1.  The  slight  difference 
in  tuning  peak  is  due  to  the  use  of  a  flat  plate  channel 
geometry  and  different  thermal  wall  boundary  conditions  in 
the  analytic  investigation. 

It  is  to  be  emphasized  that  the  correct  tuning 
condition  under  various  combinations  of  working  fluid  and 
solid  wall  material  is  a  crucial  point  in  the  design  of 
practical  devices  using  the  enhanced  thermal  pumping 
technique.  This  is  because  the  enhanced  axial  heat  flux 
which  is  usually  orders  of  magnitude  higher  than  the  axial 
heat  flow  produced  by  pure  molecular  conduction  occurs  only 
near  the  tuning  point  and  hence  at  a  fixed  radius  for  a 
fixed  value  of  w.  Once  the  tuning  point  has  been  found,  one 
can  select  the  appropriate  pipe  diameter  to  optimize  the 
heat  transfer  process. 

We  will  demonstrate  the  above  points  in  another  way. 
Fig.  4-39  shows  the  ratio  of  the  axial  heat  flux  by  pure 
molecular  heat  conduction  in  either  the  oscillating  fluid  or 
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Fig  4-3  8    Tuning  Curve  versus  Wormersley  Number 

(after  Kurzweg,   [16]  ) 
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Fig  4-39    Ratio  of  Heat  Conduction  to  Enhanced  Heat  Flux 

Versus  Wormersley  Number 
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in  the  solid  wall  to  the  enhanced  axial  heat  flux  for  both 
the  water-glass  and  mercury-steel  combinations.  These 
curves  look  like  up-side-down  tuning  curves.  They  clearly 
show  that  only  at  or  near  the  tuning  points  the  axial  heat 
flux  generated  by  the  pure  molecular  heat  conduction  either 
in  fluid  or  the  conducting  wall  is  quite  small  compared  to 
the  enhanced  heat  flux  and  hence  negligible  (about  5  lO"^ 
for  the  water-steel  wall  case,  and  about  5  10~2  for  the 
mercury-steel  wall  case) .  However,  if  one  departs  from  the 
tuning  point  this  ratio  grows  quickly  and  eventually  may  be 
greater  than  1.  The  tube  radii  corresponding  to  the  tuning 
points  (a  =  1  for  water-steel,  a  =  15  for  mercury-steel)  at 
w  =  25  rad/sec,  correspond  to  R^  =  0.02  cm  for  water-steel 
and  Ri  =  0.1  cm  for  mercury-steel . 
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CHAPTER  V 
CONCLUDING  REMARKS 


We  have  attempted  to  present  a  systematic  numerical 
investigation  of  enhanced  axial  heat  transfer  in  oscillating 
pipe  flow  for  both  insulating  wall  and  conducting  wall 
cases.  As  has  been  shown  in  Fig.  1-1,  the  pipes  which 
,  connect  both  reservoirs  are  the  most  important  elements  in 
such  enhanced  thermal  pumping  devices  and  the  performance  of 
the  oscillating  fluid  within  the  pipes  will  determine  to  a 
large  extent  if  such  devices  are  practical  for  a  given 
^  ■  application  such  as  cooling  of  micro-circuits  or  the  removal 
of  heat  from  nuclear  reactors  without  associated  convective 
mass  exchange. 

The  governing  differential  equations  derived  in 
chapter  II  are  based  on  the  assumption  that  viscous 
dissipative  terms  in  the  energy  equation  are  negligible 
compared  to  the  convective  term  UT,  and  that  the  final 
oscillating  flow  can  be  well  approximated  by  a  1-D  laminar 
oscillating  flow  with  constant  thermodynamic  properties. 

The  numerical  calculations  presented  in  this 
dissertation  obtain  solutions  by  breaking  up  the  governing 
equations  into  their  finite  difference  form  and  solving  the 
momentum  equation  with  the  Crank-Nicolson  technique,  and  the 
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energy  equations  with  a  time  dependent  Alternating  Direction 
Implicit  method  (ADI) .  a  computer  code  named  ETP  (Enhanced 
Thermal  Pumping)  was  developed  to  implement  the 
calculations,   it  is  documented  in  the  appendix. 

The  physical  mechanism  for  the  enhanced  thermal 
pumping  technique  in  oscillating  pipe  flow  is  an  interchange 
of  heat  between  the  core  slug  flow  and  the  boundary  layer 
and  the  bounding  conducting  wall.   m  fact,  the  process  acts 
very  much  as  an  enhanced  molecular  diffusion  process  in 
which  the  tidal  displacement  Ax  plays  a  role  similar  to  the 
phonon  mean  free  path  in  the  molecular  conduction  and  the 
existence  of  large  radial  temperature  gradients,  orders  of 
magnitude  higher  than  those  existing  axially,  allow  very 
large  radial   conductive  fluxes.    since  the  macroscopic 
distance  ax  is  orders  of  magnitude  larger  than  the  molecular 
mean  free  path,  it  is  not  surprising  that  axial  heat  flows 
orders  of  magnitude  larger  than  those  possible  by  conduction 
in  the  absence  of  oscillations  become  possible. 

This  enhanced  thermal  pumping  technique  is  especially 
suited  for  those  problems  where  it  is  desirable  to  transport 
large  amounts  of  heat  without  an  accompanying  convective 
mass  exchange.  The  removal  of  heat  from  radioactive  fluids 
would  appear  to  be  ideally  suited  for  a  heat  transfer  device 
based  on  this  enhanced  heat  exchange  technique. 

Three  Models  of  different  configuration  were  examined 
numerically  in  detail.   Model  1  was  employed  to  simulate 
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constant  pipe  end  temperature  conditions  as  widely  used  in 
theoretical  studies.   The  performance  of  oscillating  pipe 
flow  with  this  model  has  been  examined  numerically  and  the 
computations    well    justify    the    existing    analytical 
approximations  used.    It  is  now  clearer  that  the  tuning 
condition  and  a  large  maintained  axial  temperature  gradient 
between  pipe  ends  is  crucial  for  an  optimal  functioning 
thermal  pump.   it  was  found  that  the  enhanced  axial  heat 
flux  achievable  with  this  model  is  indeed  very  large  and  can 
exceed  by  orders  of  magnitude  the  heat  flux  possible  by  pure 
heat  conduction  between  the  pipe  ends.   The  enhanced  axial 
heat  flux  for  this  geometry  is  proportional  to  the  second 
power  of  the  tidal  displacement  for  both  the  conducting  wall 
and  insulating  wall  case  and  with  either  water  or  mercury  as 
the   working   medium   (Figs.   4-29,   and   4-30).      The 
proportionality  coefficients  between  <t>     and  ax2  have  been 
found  for  the  above  mentioned  test  cases.   Unfortunately, 
with   the   present   second   order   numerical   method,   the 
numerical  solutions  were  distorted  primarily  by  dispersion 
error  at  both  pipe  ends  where  a  discontinuity  of  the 
physical  properties  occurs  as  Wormersley  Number  and/or  tidal 
displacement  become  large. 

Model  2  was  designed  to  examine  the  properties  of  an 
alternative  thermal  pump  design  of  potential  use  in  micro- 
circuit  cooling.  such  a  configuration  was  successful  in 
transferring  large  quantities  of  heat  between  the  iso- 
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thermal  side  wall  sections  and  the  fluid  when  small  tidal 
displacements  were  used.  However,  a  limit  on  the  increase 
of  axial  heat  flux  with  increase  in  tidal  displacement 
exists  in  both  the  oscillating  and  steady  flow  cases  in  this 
model  and  the  quadratic  relationship  no  longer  holds  when 
the  tidal  displacement  becomes  large  because  of  an 
incompatibility  of  large  time-dependent  axial  heat  transfer 
in  the  thermal  pumping  process  and  the  capability  of  the 
thermal  sources  to  supply  and  remove  sufficient  heat  by 
thermal  conduction.  Nevertheless,  it  made  us  pay  more 
attention  to  a  study  of  the  heat  exchange  problem  at  the  end 
points  of  the  pipes  where  they  connect  to  fluid  reservoirs. 
Since  an  effective  enhanced  thermal  pump  produced  by 
oscillating  pipe  flow  requires  rapid  heat  supply  and  removal 
at  both  pipe  ends  any  configuration  unable  to  do  this  will 
not  allow  one  to  take  full  advantage  of  the  process. 

To  investigate  the  influence  of  the  heat  supply  and 
removal  at  the  pipe  ends.  Model  3  was  constructed  by  adding 
very  long  conductive  extensions  at  both  ends  of  the  central 
pipe  section.  This  model  simulates  the  situation  of 
oscillating  pipe  flow  for  which  no  convective  heat  exchange 
occurs  in  the  end  reservoirs  but  large  areas  for  conduction 
heat  transfer  are  made  available.  Numerical  results  show 
that  the  effective  enhanced  axial  heat  transfer  is  greatly 
weakened  over  that  achievable  by  good  convective  mixing  at 
the  pipe  ends.    The  reason  is  again  that,  even  with  an 
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enlarged  heat  source  size  (pipe  extension)  it  is  still  not 
able  to  match  the  very  fast  enhanced  axial  heat  transfer 
possible  in  the  connecting  pipe.  It  is  found  that  constant 
end  temperature  (at  x  =  5L  and  6L)  can  not  be  maintained  and 
the  quadratic  relationship  between  the  enhanced  heat  flux 
and  the  tidal  displacement  will  no  longer  hold  (Fig.  4-35) . 
However,  even  with  this  model,  the  numerical  results  show 
the  tuning  point  concept  exists.  For  water,  this  tuned 
value  occurs  at  about  Wormersley  number  a  =  1.  It  is  also 
shown  that,  under  tuned  conditions,  the  enhanced  axial  heat 
transfer  is  orders  of  magnitude  larger  than  that  possible  by 
pure  molecular  conduction.  However,  any  deviation  from  the 
tuning  condition  will  greatly  decrease  the  effectiveness  of 
the  enhanced  thermal  pumping. 

A  study  on  the  influence  of  wall  thickness  and  pipe 
diameter  on  enhanced  heat  transfer  was  performed.  The 
numerical  predictions  show  that  for  water,  small  pipe 
diameter  will  be  beneficial  and  the  optimum  wall  thickness 
should  be  about  20  percent  of  the  pipe  inner  radius  for  the 
case  considered  (i.e.,  for  inner  pipe  radius  of  R^  =  0.1cm 
the  best  outer  pipe  radius  should  be  chosen  as  R2  =  0.12cm). 
The  influence  of  the  variation  of  thermal  and  viscous 
properties  of  working  fluid  on  the  axial  heat  flux  has  also 
been  studied.  The  numerical  solutions  show  that  even  in  the 
relatively  narrow  temperature  range  from  0°C  to  100° C,  the 
axial  heat  flux  varies  more  than  150  percent. 
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Extensions  of  the  numerical  studies  considered  here 
should  include:  1)  developing  a  3-D  or  at  least  an 
axisyininetric  2-D  model  which  is  able  to  best  approximate  the 
entire  system  including  both  reservoirs  so  that  one  may  more 
accurately  study  the  effect  of  heat  exchange  within  the 
reservoirs;  2)  the  role  of  turbulence  in  such  oscillating 
flows  at  high  frequency  w  and  tidal  displacement  Ax  should 
receive  special  attention  as  well  as  the  inertia  forces  that 
are  known  to  become  large  at  higher  oscillating  frequency. 
In  fact,  a  turbulent  modelling  is  necessary  especially  if 
the  study  includes  heat  source  and  heat  sink  reservoirs;  3) 
a  consideration  of  the  variation  of  the  thermodynamic 
properties  of  working  fluid  with  temperature  and  hence 
spatial  position  should  also  be  included  in  any  extension  of 
this  research  on  enhanced  heat  transfer  by  oscillations  of 
viscous   fluids  in  pipes. 
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APPENDIX 
ETP  COMPUTER  CODE 

An  axisymmetric  code  with  designation  ETP  has  been 

developed  for  computing  the  time-dependent  Enhanced  Thermal 

Pumping  problem  with  either  oscillating  flow  (pipe  or  flat 

plate)   or  steady  flow.    The  grid  is  generated  with  a 

clustering   in   those   regions   where   high   velocity   and 

temperature  gradients  are  expected.   Both  Model  1  and  Model 

3  geometry  are  programed  in  the  current  code.    Another 

separate  code  (MIXPIP)  is  also  available  for  the  Model  2 

geometry  which  is  not  discussed  in  this  appendix. 

Description  of  Input  Variables 

The  following  statements  are  used  in  the  ETP  code: 

READ  (5,10)  NMAX,KMAX,JMAX,ICLUS,IPIPE,IPRES         '  ' 

READ  (5,10)  KMID,IREAD,JPER,ICH 

READ  (5,20)  DS,RMUK,ROU,CP,RKS,CKS 

READ  (5,20)  ALFA,PR,RL,T01,T02,AMP 

READ  (5,20)  RAD,  RADl ,RLEND,UEPS , TEPS 

READ  (5,10)  (IH(I),I=1,ICH) 

10    FORMAT  (8110) 
20    FORMAT  (8F10.5) 

where : 

NMAX  total  time  steps  within  one  period 

KMAX  total  nodal  points  on  radius 

JMAX  total  nodal  points  on  axis 
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ICLUS  =  0  . .  grid  with  clustering 

=  1  . .  uniform  Grid 
IPIPE  =  0  . .  Pipe  geometry 

=  1  . .  2-D  Flat  plate  geometry 
IPRES  =  0  . .  sinusoidal  pressure  gradient 
=  1  . .  constant  pressure  gradient 

KMID  node  index  at  interface  along  radius 

KMID  =  KMAX insulating  wall 

<  KMAX conducting  wall 

IREAD  =  0  ..  start  a  new  job 

=-2  . .  start  flow  part  only 
=  2  . .  start  heat  part  only 

(based  on  existing  flow  variables) 
=  1  . .  restart  both  flow  &  heat  field  calculation 
=  5  . .  restart  flow  field  calculation  only 

JPER number  of  period  in  the  current  run 

ICH  =  1 Model  1 

=4  ....  Model  3 

DS  distance  of  grid  next  to  wall 

RMUK  kinematic  Viscosity  of  fluid 

ROU density  of  fluid 

CP specific  heat  of  fluid 

RKS  thermal  dif fusivity  of  solid  wall 

CKS  heat  conductivity  of  wall 

ALFA  Wormersley  number 

PR  Prandtl  number 


151 

RL Length  L  of  pipe  or  flat  plate 

TOl  Temperature  in  the  hot  reservoir 

T02  Temperature  in  the  cold  reservoir 

AMP  The  amplitude  of  axial  pressure  gradient 

AMP  =  -(l./ROU)*DP/DX 

RAD inner  radius  of  pipe 

RADl  outer  radius  of  pipe 

RLEND  each  extension  portion  length 

UEPS the  allowance  velocity  residual 

TEPS the  allowance  temperature  residual 

IH(I)  node  index  of  separating  segments  in  Model  3 

Note: 

to  run  and  restart  a  job  with  this  code,  one  of  the 

procedures  has  to  be  followed: 

1)  IREAD  =  0  start  both  flow  &  heat  field 

=  1  restart  both  flow  &  heat  field 

or 

2)  IREAD  =-2  start  flow  field  only 

=2  start  heat  &  restart  flow  field 

=  1  restart  both  flow  &  heat  fields 

A  Sample  Input  Data  and  Batch  Command  File 

The   following   sample   input   was   used   in   the 

computation  with  Model  3  (water-glass) 

1001        15       101  0         0         0 

15          0          2  4 

0.40000  0.01006    0.99820  4.18181  0.20260    0.73000 

10.00000  7.02000   20.00000  1.00000  0.00000  29270.600 

0.10000  0.15000  100.00000  0.00500  0.00500 

1         21         81  101 
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The  following  Command  file  was  used  to  start  the  run 

of  a  new  case: 

$  ASSIGN  SAMPLO.INP  FOR005 

$  ASSIGN  SAMPLO.DAT  FOR006 

$  ASSIGN  FILE7.DAT  FOR007 

$  ASSIGN  FILE8.DAT  FOR008 

$  ASSIGN  FILE11.DAT  FOROll 

$  R  EPT 

To  restart  the  previous  run,  the  following  command 

file  was  used: 

$  RENAME  FILE7.DAT  FILE9.DAT 

$  RENAME  FILE8.DAT  FILE10.DAT 

$  RENAME  FILE11.DAT  FILE12.DAT 

$  ASSIGN  SAMPLl.INP  FOR005 

$  ASSIGN  SAMPL1.DAT  FOR006 

$  ASSIGN  FILE7.DAT  FOR007 

$  ASSIGN  FILE8.DAT  FOR008 

$  ASSIGN  FILE9.DAT  FOR009 

$  ASSIGN  FILE10.DAT  FOROIO 

$  ASSIGN  FILE11.DAT  FOROll 

$  ASSIGN  FILE12.DAT  FOR012 

$  R  EPT 

Note,  the  name  of  the  data  files  could  be  changed  according 

user's  taste. 

Function  of  Subroutine  in  ETP  Code 

In  the  ETP  code,  the  Program  MAIN  controls  the  main 
loop  in  the  computation,  including  the  input  and  output  data 
control,  restart  control.  Also  it  sets  up  the  initial 
condition  according  to  the  purpose  of  investigation.  The 
ETP  code  includes  22  subroutines.  The  function  of  each 
subroutine  is  listed  below: 
Subroutine  COMMENT 

This  subroutine  gives  the  description  of  the  input 
and  output  variables  listed  in  previous  section. 


^H  '^-  '        153  ,:   '•    ■     .■   ^-  -' 

Subroutine  GRID        ^ 

f'--""   ■■■'  * 

f"l--'-  "  This  subroutine  is  used  to  generate  either  a  uniform 

or  a  clustered  grid  system  in  the  physical  domain  for  Model 
1  or  Model  3;  it  calls  the  subroutine  CLUS  to  cluster  the 
grid  lines  in  the  vicinity  of  wall  and  the  central  part  of 

•J;      pipe  (Model  3)  .   It  also  calls  subroutine  METPIP,  METWAL  to 
compute  the  derivatives  terms  in  the  transformation. 

;\      Subroutine  CLUS  (ICLUS.K2) 

This  subroutine  is  called  by  GRID  to  compute  the 

';^  ,   clustering  grid  distribution  with  equation  (3-52). 
Subroutine  METPIP 

This  subroutine  is  called  to  compute  the  derivatives 
of  the  coordinate  transformation  in  the  inner  pipe  area. 

Such  as  dx/d^,    dx/dri     used  in  equations  from  (3-6)  to 

(3-16)  as  well  as  in  the  boundary  condition. 
Subroutine  METWAL 

Similarly,  subroutine  METWAL  is  called  with  the  same 
equations  to  compute  the  derivatives  in  transformation  in 
the  solid  wall,  if  a  conducting  wall  is  assumed. 

:      Subroutine  FLOW  fN) 

This  subroutine  is  called  to  compute  the  velocity  in 
the  flow  field  at  each  time  step  with  the  Crank-Nicolson 
Method.   Equation  (3-20)  or  the  matrix  form  (3-28)  will  be 
solved  in  this  subroutine. 
Subroutine  TIDAL 

1^5  This  subroutine  is  called  to  Compute  the  Lagrangian 
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displacement  (2-33),  tidal  displacement (2-32)  and  the  phase 
lags  which  are  relative  to  the  exciting  pressure  gradient 
phase. 
Subroutine  COFLOW 

This  subroutine  is  called  to  compute  the  coefficients 
in  the  momentum  equation  at  each  time  step,  i.e.,  computing 
equations  (3-7)  and  (3-8) . 
Subroutine  RHSFLO  (IJ) 

The  subroutine  RHSFLO  (N)  is  called  to  compute  the 
update  right  hand  side  terms  in  the  momentum  equation  at 
each  time  step,  i.e.,  equations  from  (3-21)  to  (3-25).  = 
Subroutine  TRIP 

A  subroutine  TRID  (borrowed  from  [2])  is  to  solve  the 
tridiagonal  system  algebra  equations,  it  is  widely  used  to 
compute  the  update  velocity  and  temperature. 
Subroutine  EPSIL 

Function  EPSIL  (borrowed  from  code  "GRIDGEN")  to  find 
the  €  value  used  in  equation  (3-52)  with  Newton-Raphson 
root-finding  technique. 
Subroutine  ADISU3 

This  subroutine  is  called  using  ADI  method  to  compute 
the  temperature  distribution  in  oscillating  pipe  flow  in 
Model  3.  The  insulating  wall  boundary  condition  is  assumed 
in  the  central  pipe  (5L  <  x  <  6L)  with  fixed  temperature  at 
the  far  ends  of  both  extension  pipe.  Equations  from  (3-29) , 
(3-35) ,  or  the  matrix  from  (3-33) ,  (3-39)  are  solved. 
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Subroutine  COTPIP 

This  subroutine  is  called  by  ADISU3  to  calculate  the 
coefficients  in  the  energy  equation  of  pipe  flow  (3-41)  at 
each  time  step,   which  are  coupled  by  the  update  velocity. 
Subroutine  RHSINT  (N. II. LINE) 

This  subroutine  is  called  by  ADISU3  to  compute  the 
update  right  hand  side  terms  in  the  energy  equation  of  pipe 
flow  with  the  insulating  wall  case,  namely,  solving  (3-31) , 
(3-37) . 
Subroutine  FLUXT  fNN) 

Subroutine   FLUXT   (NN)   is   called   to   compute  the 
dimensionless  enhanced  axial  heat  flux  as  well  as  the  heat 
flux  by  pure  molecular  conduction  either  in  the  wall  or  the 
pipe  fluid,  namely,  equations  (2-34)  to  (2-37). 
Subroutine  OUT  fNTID) 

The  subroutine  OUT  (NTID)  is  called  to  compute  and 
output  the  dimensional  heat  flux,  tidal  displacement,  and 
the  average  temperature  (dimensionless) . 
Subroutine  ADIC03 

Subroutine  ADIC03  is  similar  to  ADISU3,  but  for  the 
conducting  wall  case. 
Subroutine  COTWAL 

This  subroutine  is  called  to  calculate  the 
coefficients  in  the  wall  conduction  equation  at  each  time 
step  by  using  equation  (3-42) . 
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Subroutine  RHSCOT  TN. II. LINE) 

Subroutine  RHSCOT  is  similar  to  RHSINT,  but  for  the 
conducting  wall  case. 
Subroutine  RHSWAL  rN.II.LINE^ 

Subroutine  RHSWAL  (N, II, LINE)  is  called  to  compute 
the  right-hand-side  terms  in  the  wall  conduction  equation 
with  (3-32)  and  (3-38) . 
Subroutine  ADISUl  (N)  '  ' 

Subroutine  ADISUl  is  similar  to  ADISU3 ,  but  for  the 
case  of  Model  1. 
Subroutine  ADICOl  CN) 

Subroutine  ADICOl  is  similar  to  ADISUl,  but  for  the 
conducting  wall. 


REFERENCES  -^ 

1.  Abramowitz,  M. ,  and  Stegun,  I.  A. (editors) ,  "Handbook  of 
Mathematics  Functions  with  Formulas,  Graphys,  and  Mathe- 
matical Tables",  Dover  Publications,  Inc.,  New  York 

(1972)  .  v. 

2.  Ames,  W.  F. ,  "Numerical  Method  For  Partial  Differential 
Equations,"  Academic  Press,  INC. 

3.  Anderson   D.  A.,  Tannehill,  J.  C. ,  and  Fletcher,  R.  H. , 

'Computational  Fluid  Mechanics  and  Heat  Transfer  " 

McGraw-Hill,  New  York  (1984).  '  '    r*-'\ 

^'    mf^^'  ?•'  '"^^^  Dispersion  of  Solute  in  Pulsating  Flow 
Through  a  Tube,"  Proceedings  of  Royal  Society,  London, 
A  259,  pp. 370-376  (1960).  ' 

5.  Arpaci  V.,  and  Larsen  P.,  "Convective  Heat  Transfer  " 
Prentice  Hall,  Englewood  Cliffs,  NJ  (1984). 

6.  Baldwin  B.  S.,  and  Lomax,  H. ,  "Thin  Layer  Approximation 
and  Algebraic  Model  for  Separated  Turbulent  Flow  "  AIAA 
16th  Aerospace  Sciences  Meeting  (Jan.,  1978)  pp. 78-257 

7.  Bohn   D.  J.,  Miyasaka,  K. ,  Marchak,  B.  E. ,  Thompson,  W. 
K.,  Froese,  A.  B. ,  and  Bryan,  A.  C. ,  "Ventilation  by 
High-Frequency  Oscillation,"  Journal  of  Appl.  Phvs 
Vol.  48,  pp. 710-716  (1980).  ^   ' 

8.  Bowden,  K.  F. ,  "Horizontal  Mixing  in  the  Sea  Due  to 
Shearing  Current,"  Journal  of  Fluid  Mech. ,  Vol.  21 
p. 84  (1965).  ' 

9.  Chapman,  A.  J.,  "Heat  Transfer,"  Macmillan  Publishing 
Co.,  New  York  (1984).  iJJ-J.t.uing 


10, 


Chatwin,  P.  c,  "On  the  Longitudinal  Dispersion  of 
Passive  Contaminant  in  Oscillatory  Flow  in  Tubes" 

^975^  °^  ^^""^"^  ^^''^''    ^°^'  ^^'  ^^""^  ^'  PP-513-527 

"'■■'■*  «^!^®^^'  ?'  ^'    ^•'    "Analysis  of  Heat  and  Mass  Transfer  " 
McGraw-Hill,  New  York  (1972).  ' 


-^J 


157  .     .^ 


/  . 


,  158        -  .  -'   ■■'  -   r  ■    ■■;^-; 

12.  Jaeger,  M.  J.,  and  Kurzweg,  U.  H. ,  "Determination  of  the, 
Longitudinal  Dispersion  Coefficient  in  Flows  Subjected 

to  High-Frequency  Oscillations,"  Phys.  Fluid,  Vol.  26        ; 
No.  6,  pp. 1380-1382  (1983). 

13.  Joshi,  C.  H.,  Kamm,  R.  D. ,  Drazen,  J.  M. ,  and  Slutsley, 
A.  S.,  "An  Experimental  Study  of  Gas  Exchange  in  Laminar 
Oscillating  Flows,"  Journal  of  Fluid  Mech.,  Vol.  133, 
pp. 245-254  (1983). 

14.  Kaviany,  M. ,   "Some  Aspects  of  Heat  Diffusion  in  Fluids 
by  Oscillation,"  International  Journal  Heat  Mass 
Transfer,  Vol.  29,  pp. 2002-2006  (1986). 

15.  Kurzweg,  U.  H. ,  "Enhanced  Heat  Conduction  in  Fluids 
Subjected  to  Sinusoidal  Oscillations,"  Journal  of  Heat 
Transfer,  Vol.  107,  pp. 459-462  (1985). 

16.  Kurzweg,  U.  H. ,  "Enhanced  Heat  Conduction  in  Oscillating 
Viscous  Flows  Within  Parallel-Plate  Channels,"  Journal 
of  Fluid  Mech.,  Vol.  156,  pp. 291-300  (1985). 

17.  Kurzweg,  U.  H. ,  "Temporal  and  Spatial  Distribution  of 
Heat  Flux  in  Oscillating  Flow  Subjected  to  an  Axial 
Temperature  Gradient,"  International  Journal  Heat  Mass 
Transfer,  Vol.  29,  No.  12,  pp. 1969-1977  (1986). 

18.  Kurzweg,  U.  H. ,  "Enhanced  Diffusional  Separation  in 
Liquids  by  Sinusoidal  Oscillations,"  to  appear  in 
Separation  Sciences  and  Technology,  in  press. 

19.  Kurzweg,  U.  H. ,  and  Jaeger,  M.  J.,  "Diffusional 
Separation  of  Gases  by  Sinusoidal  Oscillations,"  Phys. 
Fluids,  Vol.  30,  pp. 1023-1026  (1987). 

20.  Kurzweg,  U.  H. ,  and  Jaeger,  M.  J.,  "Tuning  Effect  in  Gas 
Dispersion  under  Oscillatory  Conditions,"  Phys.  Fluids 
Vol.  29,  pp. 1324-1326  (1986). 

21.  Kurzweg,  U.  H.  and  Lindgren,  E.  R.  "Enhanced  Heat 
Conduction  By  Oscillatory  Motion  of  Fluids  in  Conduits," 
A  Research  Proposal  to  the  Fluid  Dynamics  and  Heat 
Transfer  Program,  The  National  Science  Foundation, 
under  Contract  Number  CBT-8611254  (1986). 

22.  Kurzweg,  U.  H. ,  and  Zhao,  L.  D. ,  "Heat  Transfer  by  High- 
Frequency  Oscillations:  A  New  Hydrodynamic  Technique 
for  Achieving  Large  Effective  Thermal  Conductivities," 
Phys.  Fluids,  Vol.  27,  pp. 2624-2627  (1984).. 


■  -A 


,H  ■ 


-W"' 
^.''M^ 

^^. 


159 

23.  Li,  C.  P.,  "A  Finite  Difference  Method  for  Solving 

Unsteady  Viscous  Flow  Problems,"  AIAA  Journal,  Vol.  23, 
No. 5  (1985). 

24.  MacCormack,  R.  W. ,  "Current  Status  of  Numerical 
Solutions  of  the  Navier-Stokes  Equations,"  AIAA  23rd 
Aerospace  Sciences  Meeting,  Reno,  Nevada,  AIAA-85-0032 
Jan.  (1985) 

25.  Mastin,  C,  W. ,  "Error  Induced  by  Coordinate  System," 
in  J.  F.  Thomson  (editor) ,  Numerical  Grid  Generation, 
North-Holland,  New  York,  pp. 31-40  (1982). 

26.  Merkli,  P.,  and  Thomann,  H. ,  "Transition  to  Turbulence 
in  Oscillating  Pipe  Flow,"  Journal  of  Fluid  Mech. ,  Vol. 
68,  pp. 567-575  (1975). 

27.  Ommi,  M.,  Iguchi,  M. ,  and  Urahata,  I.,  "Flow  Patterns 
and  Frictional  Losses  in  An  Oscillating  Pipe  Flow," 
Bulletin  of  the  JSME,  Vol.  25,  No.  202.  April  (1982). 

28.  Pulliam,  T.  H. ,  "Efficient  Solution  Methods  for  the 
Navier-Stokes  Equations"  Lecture  Notes  for  the  Von 
Karman  Institute  (Jan.,  1986). 

29.  Schlichting,  H. ,  "Boundary  Layer  Theory"   McGraw-Hill, 
New  York  (1972) . 

30.  Shah,  V.  L. ,  "COMMIX-IB:  A  Three-Dimensional  Transient 
Single-Phase  Computer  Program  for  Thermal  Hydraulic 
Analysis  of  Single  and  Multicomponent  Systems,"  Vol.  1: 
Equations  and  Numerics,  Available  from  Superintendent  of 
Documents  U.  S.  Government  Printing  Office,  P.  0.  Box 
37082,  Washington,  D. C. 20013-7982 ,  Sep.,  (1985). 

31.  Shih,  T.  M.,  "Numerical  Heat  Transfer,"  Hemisphere 
Publishing  Corp.,  Washington  (1984). 

32.  Taylor,  F.  R.  s.,  "Dispersion  of  Soluble  Matter  in 
Solvent  Flowing  Slowly  Through  a  Tube,"  Proc.  Roy.  Soc, 
pp. 186-203  (1953). 

33.  Thompson,  J.  F. ,  "Grid  Generation  Technique  in 
Computational  Fluid  Dynamics",  AIAA  Journal,  Vol.  22 
Nov. 11,  pp. 1506-1519,  (1984). 

34.  Thompson,  J.  F. ,  Thames,  F.  C.  and  Mastin,  C.  W. , 
"Automatic  Numerical  Generation  of  Body-Fitted 
Curvilinear  Coordinate  System  for  Field  Containing  Any 
Number  of  Arbitrary  Two-Dimensional  Bodies,"  Journal  of 
Computational  Physics.,  Vol.  15,  pp. 299-319  (1974). 


160 

35.  Thompson,  J.  F. ,  Warsi,  Z.  U.  A.,  and  Mastin,  C.  W. , 
"Boundary-Fitted  Coordinate  Systems  for  Numerical 

Solution  of  Partial  Differential  Equations A  Review", 

Journal  of  Computational  Physics.,  Vol  47.,  pp. 1-108 
(1982)  . 

36.  Uchida,  S.,  "The  Pulsating  Viscous  Flow  Superposed  on 
the  Steady  Laminar  Motion  of  Incompressible  Fluid  in  a 
Circular  Pipe,"  ZAMP,  Vol.  VII,  pp. 402-421  (1956). 

37.  Warming,  R.  F. ,  and  Beam,  R.  M. ,  "On  the  Construction 
and  Application  of  Implicit  Factored  Schemes  for 
Conservation  Laws,"  SIAM-AMS  Proceedings,  Vol  11, 

pp. 85-129  (1978) . 

38.  Watson,  E.  J.,  "Diffusion  in  the  Oscillation  Pipe  Flow," 
Journal  of  Fluid  Mech.,  Vol.  133,  pp. 233-244  (1983). 


BIOGRAPHICAL  SKETCH 

Guo-Jie  Zhang  was  born  on  June  26,  1941,  in  Beijing, 
China,  and  spent  his  childhood  and  youth  in  that  beautiful 
capital  city.    After  graduating  from  Beijing  2lth  nigh 
School  in  1960  he  entered  the  University  of  Science  and 
Technology  of  China,  Beijing,  China.    After  5  years  of 
college   education   in   the   Department   of  Mechanics,   he 
graduated  in  July,  1965.    since  then  he,  as  an  aircraft 
structure  strength  analysis  engineer,  has  been  engaged  at 
the  Aircraft  Research  and  Development  Institute,  Nanchang 
Manufacturing  Aircraft  Company,  China.   until  the  spring  of 
1981.   During  that  period  he  participated  in  the  research  on 
and   development   of   several   aircraft   with   his   main 
contributions  being  in  the  area  of  the  studies  with  finite 
element  methods.   m  May,  1981,  he  was  selected  to  study 
abroad  by  the  Ministry  of  Aviation  Industry  of  China  and 
after  one  year's  English   language  training  at  Nanjing 
Aeronautical   Institute,   he   was   sent   to   study   in   the 
Department   of   Aerospace   Engineering,   Mechanics,   and 
Engineering  Science,  University  of  Florida  as  a  visiting 
scholar,  and  later,  in  May,  1983,  transferred  to  graduate 
student  status.    He  received  his  Master  of  Engineering 
degree  in  April,  1986,  and  expects  to  receive  the  degree  of 

161 


162 
Doctor  of  Philosophy  in  the  spring  of  1988.   He  is  married 
to  Chun-hua  Shi  and  has  two  children,  Yan  Zhang  and  Wei 
Zhang.   Mr.  Guo-jie  Zhang  is  a  member  of  the  Aeronautics 
Society  of  China. 


•^ 


>  -  ^^    i 


I  certify  that  I  have  read  this  study  and  that  in  my 
opinion  it  conforms  to  acceptable  standards  of  scholarly 
presentation  and  is  fully  adequate,  in  scope  and  quality,  as 
a  dissertation  for  the  degree  of  Doctor  of  Philosophy. 


-_■  ^  ' 


^/h'cA  ?/'    ^k2^a 


Ulrich  H.  Kurzweg,  Chairman 
Professor  of  Aerospace  Engineering, 
Mechanics,  and  Engineering  Science 


I  certify  that  I  have  read  this  study  and  that  in  my 
opinion  it  conforms  to  acceptable  standards  of  scholarly 
presentation  and  is  fully  adequate,  in  scope  and  quality,  as 
a  dissertation  for  the  degree  of  Doctor  of  Philosophy. 


^ 


E,  Rune  Lindgren 

Professor  of  Aerospace  Engineering, 

Mechanics,  and  Engineering  Science 


I  certify  that  I  have  read  this  study  and  that  in  my 
opinion  it  conforms  to  acceptable  standards  of  scholarly 
presentation  and  is  fully  adequate,  in  scope  and  quality,  as 
a  dissertation  for  the  degree  of  Doctor  of  Philosophy. 


Lawrence  E.  Malvern 

Professor  of  Aerospace  Engineering, 

Mechanics,  and  Engineering  Science 


J^i 


I  certify  that  I  have  read  this  study  and  that  in  my 
opinion  it  conforms  to  acceptable  standards  of  scholarly 
presentation  and  is  fully  adequate,  in  scope  and  quality,  as 
a  dissertation  for  the  degree  of  Doctor  of  Philosophy. 


f 


"^=1 


vJ  i'-nj^ 


David  W.  Mikolaitis 
Assistant  Professor  of  Aerospace 
Engineering,  Mechanics,  and 
Engineering  Science 


I  certify  that  I  have  read  this  study  and  that  in  my 
opinion  it  conforms  to  acceptable  standards  of  scholarly  "* 
presentation  and  is  fully  adequate,  in  scope  and  quality,  as  •• 
a  dissertation  for  the  degree  of  Doctor  of  Philosophy. 


TU.V   \VWv.a-^^  VCXy^^ 


Arun  K.  Varma 

Professor  of  Mathematics 


This  dissertation  was  submitted  to  the  Graduate  Faculty  of 
the  College  of  Engineering  and  to  the  Graduate  School  and 
was  accepted  as  partial  fulfillment  of  the  requirements  for 
the  degree  of  Doctor  of  Philosophy. 

April  1988 


.MM^ 


Dean,  College  of  Engineering 


Dean,  Graduate  School 


UNIVERSITY  OF  FLORIDA 


3  1262  08556  7963 


